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PREFACE 


This  is  the  second  of  two  volumes  reporting  on  work  performed 
under  contract  to  the  Defense  Advanced  Research  Projects  Agency  from 
October  1975  through  November  1978.  This  volume  contains  supporting 
material  for  Vol.  I,  and  is  organized  into  six  separate  appendices. 
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APPENDIX  A 


STRESS- STRAIN  RELATIONS 


The  constitutive  equations  developed  for  this  program 

have  been  cast  in  a form  that  currently  enjoys  widespread  use 

A1 

for  geologic  materials.  The  calculation  of  mean  stress,  P, 
and  deviatoric  stress,  , are  separated.  These  calculations  are 
further  subdivided  to  account  for  the  presence  (or  absence)  of 
shear  failure.  In  the  first  part  of  this  appendix,  attention 
is  confined  to  just  one  portion  of  the  complete  constitutive 
model,  namely,  calculation  of  the  mean  stress  in  the 
absence  of  shear  failure.  In  that  case  P is  computed  from 
formulas  that  collectively  define  a "hydrostat"  — the  set  of 
all  states  that  can  be  reached  by  purely  hydrostatic  loading 
and  unloading.  Calculation  of  the  deviatoric  stress,  which 
together  with  P fully  defines  the  stress  tensor,  the  develop- 
ment of  shear  failure,  and  the  calculation  of  P when  shear 
failure  occurs  are  described  in  later  sections  of  this  appendix. 


Al, 


Mean  Stress  in  the  Absence  of  Shear  Failure 


In  the  absence  of  shear  failure,  the  mean  stress  P is 

written  as  the  sum  of  a contribution,  P , due  to  vapor,  and  a 

o 

contribution,  P , from  condensed  material: 
s 


P + P 
g s 


(Al) 


In  ATI's  constitutive  equations  for  hard  rock,  the  dependence 

of  P on  density  and  specific  internal  energy  is  expressed  in 
8 A2 

the  well-known  Tillotson  form  (Section  A1.2).  The  calculated 

of  P is  not  so  simple.  The  method  of  constructing  a hydrostat 
s 

for  a given  geologic  material  now  used  routinely  at  ATI  is  out- 


A1 


mmm 


lined  in  Section  Al.l.  Thermal  expansion  of  the  solid,  not 

routinely  accounted  for  in  ATI's  procedure  for  calculating  P , 

s 

is  taken  into  account  in  the  work  reported  here. 

Al.l  The  Variation  of  Mean  Stress  with  Excess  Compression 
in  Cold  Solid  (No  Shear  Failure) 

On  hydrostatic  loading  and  subsequent  unloading,  geologic 

solids  generally  experience  irreversible  volume  changes  — a hys- 

teretic  process  termed  "compaction".  In  all  but  the  hardest 

rocks,  compaction  leads  to  rapid  attenuation  of  stress  waves, 

kinetic  energy  is  dissipated  and  locally  heats  the  rock.  It 

is  important  to  model  such  behavior  faithfully  in  calculating 

P , the  contribution  made  to  the  mean  stress  by  the  condensed 
s 

component  of  a geologic  medium.  To  that  end,  it  is  useful  to 

define  an  "excess  compression",  , which  is  equivalent  to  the 

volume  strain  of  a material  element;  specifically,  p,  is  equal 

to  p/p  -1,  where  p is  the  density  of  a given  material  in  a 
o o 

reference  state  (zero  stress)  and  p is  its  density  in  an  arbi- 
trary  state  of  stress.  By  allowing  P to  vary  with  y,  different- 
ly  in  loading  than  in  unloading,  as  will  now  be  discussed  at 
some  length,  permanent  decreases  in  volume  can  be  provided  for 
in  the  constitutive  model. 

4 

If  the  mean  stress  Pg  is  increased  continuously,  material 
traverses  a set  of  "initial  loading"  or  "virgin  loading"  states 
represented  as  a single  curve  in  the  P-y  plane.  Experimental 
data  show  that  a material  containing  cracks  and/or  pores  is  con- 
siderably  more  compactible  (has  a lower  bulk  modulus)  at  low 
values  of  mean  stress  than  it  has  at  high  mean-stress  values, 

Typically,  as  pores  and  open  cracks  close  the  bulk  modulus,  K, 
increases.  Details  of  the  structure  will  affect  the  variation 
of  the  bulk  modulus  with  pressure.  In  some  cemented  materials 
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such  as  sandstone  the  bulk  modulus  may  initially  increase  with 
pressure,  decrease  as  the  cementation  breaks  down,  and  then 
increase  during  subsequent  pore  closure.  In  any  event,  however, 
intergranular  spaces  gradually  disappear  and  K approaches  its 
original  value.  On  further  compression,  the  bulk  modulus  ra- 
pidly increases,  asymptotically  approaching  some  intrinsic  value, 

K , characteristic  of  the  consolidated  material 
m 

If  the  mean  stress  is  increased  sufficiently  and  then  de- 
creased, one  finds  experimentally  that,  as  unloading  begins,  the 
bulk  modulus  has  a value  greater  than  or  equal  to  that  last  en- 
countered in  loading;  on  further  unloading,  the  effective  bulk 
modulus  monotonically  decreases  from  its  initial  unloading  value 
to  some  lower  value  at  P=0.  The  average  slope  of  a given  unload- 
ing path  in  the  P-n  plane,  however,  is  generally  greater  than 
that  of  the  loading  path  that  precedes  it,  as  a result  of  a net 
decrease  in  volume  in  the  load-unload  process.  If  the  material 
is  reloaded  after  unloading,  then,  to  a satisfactory  approxima- 
tion, the  unloading  path  is  retraced.  Accordingly,  the  process 
of  unloading  and  reloading  is  considered  elastic  up  to  the  lar- 
gest value  of  mean  stress  reached  on  initial  loading. 

Since  the  peak  stress  reached  on  initial  loading  is  arbi- 
trary, the  process  of  unloading  and  reloading  is  described  not 

by  a single  P path,  but  by  a family  of  P -n  curves.  Experi- 
s s 

mental  data  indicate  that  in  the  P-u  plane,  the  unloading/re- 
loading paths  of  most  materials  fall  into  two  classes.  For  the 
simpler  class,  the  shapes  of  the  unloading  paths  are  independent 
of  the  maximum  excess  compression,  umax>  achieved  by  a material 
element;  the  shapes  of  the  unloading  paths  depend  in  the  more 
general  case  on  the  maximum  compression  to  which  the  element  of 
material  has  been  subjected. 
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Figure  A1  presents  a qualitiative  picture  of  the  kind  of 
P-n  relation  implied  by  the  general  model  just  described.  A more 
precise  formulation  is  as  follows:  if  ii*u  for  an  element  of 
material,  then  the  path  of  initial  loading  (also  called  the 
"virgin  loading  path")  is  used  to  determine  the  mean  stress,  and 

then  the  element  is  either 


ti  is  updated, 
max  r 


However,  if  u<n  , 
’ max* 


being  unloaded  or  reloaded,  and  the  value  of  nmax  uniquely  deter- 


mines the  P -n  path  which  the  element  must  follow 

b 


Experimental 

data  show  that  until  all  cracks  and  voids  are  closed,  the  per- 
manent compaction,  nz  (the  value  of  u at  zero  mean-stress  on  an 
unloading  path),  increases  with  the  peak  value  attained  by  the 
mean  stress,  i.e.,  the  material  experiences  a continuous  and  ir- 
reversible decrease  in  porosity  as  the  mean  stress  is  increased. 
In  practical  terms,  however,  nz  attains  a finite  limiting  value, 


^zm’ 


and  y,  does  not  change  thereafter;  the  unload-reload  path 
z 


traced  in  the  P-y,  plane  on  unloading  from  Pm  is  designated  the 
"limiting  unloading  path". 

For  numerical  convenience,  y is  defined  as  a function  of 

* z 

y , while  P is  most  easily  expressed  in  terms  of  y and  y ; 
max  z 

thus  P can  be  determined  along  an  unload-reload  path  knowing 
y and 
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Thermal  Effects:  Expansion  and  Vaporization  of  Condensed 
Material 


For  the  media  treated  during  this  program,  the  contribution 
of  vapor  to  the  mean  stress,  P , was  computed  according  to  the 

O 

following  formula: 

0 ; e<e 
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al  + 


f2  + =i=Tl)  “PM 
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(A2) 


o(e-ev) 


e*e 
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VIRGIN  LOADING  PATH 


~ value  of  u at  which  an  arbitrary  unloading  path 
intersect!  the  virgin  loading  path 


li  ~ Umax  f°r  the  limiting  unloading  path;  also,  value 

I®  of  u above  which  no  further  permanent  compaction 

is  allowed 

U ~ permanent  set  for  an  arbitrary  unloading  path 

I o 

u ~ u0  f°r  the  limiting  unloading  path 


Figure  Al.  Typical  Behavior  of  a Compactlble  Material 
Under  Hydrostatic  Loads 


where  T1  and  e are  the  compression,  p/pq,  and  specific  internal 
energy  of  the  material;  er  is  equal  to  five  for  expanded  material 
(T)<1)  and  zero  when  the  material  is  compressed  (Till).  The 
characteristic  energy,  ev>  will  be  discussed  later  in  this  sec- 
tion. The  values  of  a^,  a2>  b,  etc.,  were  chosen  for  the  fol- 
lowing reasons : 

(a)  In  the  limit  of  low  density,  Eq.  (A2)  reduces  to  the 
equation  of  state  of  an  ideal  gas  with  internal  degrees  of  free- 
dom; a^  assumes  the  role  of  the  quantity  y-1,  where  y is  the 
familiar  heat  capacity  ratio.  Although  it  is  known  that  y-1 

varies  with  P,  p,  and  e,  a constant,  representative  value  of  a. 

A3  1 

has  been  adopted  here 

(b)  The  term  proportional  to  b in  Eq.  (A2)  becomes  negli- 
gible when  e greatly  exceeds  eQ  (eQ  is  identified  roughly  with 
the  material's  vaporization  energy).  Since  e is  also  much 
greater  than  ev  and  ■!■!■«  1 , the  right-hand  member  of  Eq.  (A2)  ap- 
proaches (a^  + a2>pe  and  (a.  + a^)  was  set  equal  to  .5  in  accord 

A 1 A ^ A A 

with  previous  experience  ’ * 


(c)  The  maximum  value  of  the  variable  ev  should  be  roughly 
equal  to  the  melting  energy  of  hard  rock.  Melting  energy  is 
known  to  vary  with  pressure.  Based  on  melting  energy  data  for 
materials  similar  to  those  considered  here,  the  variable  e is 
defined  as 


where 


MAX(0 . 0 , MIN  (e  , e + n)) 

' ’ v mm’  mo  dii 


et(l  + f) 


C1  * 1 VT 

Ti 


- : -v  ev  IH*- 


e is  the  maximum  melting  energy,  K,  is  the  bulk  modulus  on 
mm  i 

the  limiting  unloading  curve  at  Ps“0,  and  Cv  is  the  specific 
heat  (assumed  constant). 

To  account  for  thermal  expansion  of  condensed  material, 
a simple  approximate  procedure  was  used.  A modified  excess 
compression,  n',  was  defined  as  follows: 

n'  = H + 3EMin(e,ev)]  (A7) 

where  the  constant  B was  assigned  the  value  ®yCv  on  the  basis 
of  thermodynamic  data  for  both  the  coefficient  of  volumetric 
thermal  expansion  («v)  and  the  specific  heat  (Cv) . Obviously, 
heating  of  solid  material  causes  to  exceed  n.  If  the  mean 
stress  P is  a known  function  of  n(say  fQ(u))  in  states  for  which 
e * 0,  then  accounting  for  heating  by  substituting  n'  for  ^ in 
fQ(u) » has  the  same  effect  on  P as  an  increase  in  excess  com- 
pression at  e = 0.  Furthermore,  the  approximation  P(n ,e)*fQ(n ' ) 
is  accurate  to  the  extent  that  P varies  linearly  with  p,  and  e 
over  a range  of  states  about  standard  conditions.  This  procedure 
was  followed  for  each  material  treated  here:  \i  was  replaced  by  p,' 
in  the  function  fit  to  experimental  data  defining  the  quasistatic 
variation  of  mean-stress  with  excess  compression.  Heating  of  the 
material  at  any  given  density  then  produces  greater  compressive 
values  of  mean  stress  than  would  be  found  in  cold  solid  at  that 


density 


Al.3  Mean  Stress  on  Unload/Reload  Paths 


Curvature  of  the  P-u  paths  for  unload/reload  conditions 


by  adopting  a hydrostat  of  the  functional  form 


The  derivative  with  respect  to  x is  given  by 


In  general,  the  argument  x on  the  unload-reload  curve  is 
given  by 


where  u'  is  given  by  Eq.  (A7)  and  nz  represents  the  value  that 

would  have  if  the  material  were  unloaded  from  P = P (u'  )«P 

s v>max  max 

to  P ■ 0.  and  ^ are  bulk  moduli  on  'large'  and  'small'  ex- 
cess compression.  The  parameter  x*  controls  the  rate  of  transi- 
tion from  the  bulk  modulus  K.  to  K0. 


A simple  model  of  hydrostatic  tensile  behavior  is  used. 

Linear  behavior  is  assumed  for  P <P  <0.  However,  when  P <P.  , 

t s ’ st 

the  pressure,  bulk  modulus  and  shear  modulus  are  all  set  to 
zero.  This  is  often  referred  to  as  a 'Tension  Cut-off'  model. 
Although  these  functional  forms  provide  a general  description 
of  the  hydrostat  for  unloading  paths  and  virgin  loading  paths 
for  many  media,  description  of  hydrostatic  behavior  for  any  one 
material  may  require  some  adjustment  to  be  made.  The  specific 
equations  used  for  each  of  the  four  media  treated  here  are  given 


in  Sections  (A4)  through  (A7). 


, 
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Deviatoric  Stress  and  Shear  Failure  Relations 


A complete  description  of  the  mechanical  behavior  of 


solids  requires  a specification  of  the  deviatoric  stress,  • 


In  any  such  specification  it  is  necessary  to  provide  for  both 
elastic  and  inelastic  deformation  of  material.  Equations  for 
that  purpose  are  presented  below.  Apart  from  changes  in  the 
values  of  constants,  those  equations  were  used  throughout  the 
present  program  to  compute  increments  in  deviatoric  stress . 


A2.1  Elastic  Deviatoric  Stress  Calculation 


Unless  shear  failure  occurs,  the  incremental  deviatoric 

A5 

stress  satisfies  the  following  relation 


da! 


ij 


-2G  de! 


ij 


where  da^  and  de!^  are  increments  in  the  deviatoric  stress 


and  strain,  respectively,  and  the  shear  modulus,  G,  decreases 
with  increasing  internal  energy  e according  to  the  formula. 


where  G(n  ) is  a function  which  varies  from  material  to  ma- 
max 


terial  as  defined  in  Sections  (A4)  through  (A7),  em  is 
9G  1 

■jY  * TF^r,  and  e«>  melting  energy,  is  defined  in  Eq.  (A3). 


pc;° 


A2.2  Shear  Failure 


(A13) 


(A14) 


During  plastic  deformation,  i.e.,  in  states  of  shear  fail- 


ure, the  calculation  of  the  incremental  stress  deviator,  da|^ , 


and  of  the  mean  stress  increment  dP,  are  based  on  the  method  of 

A5 

the  plastic  potential  . A general  description  of  the  methods 


used  to  obtain  stresses  during  shear  failure  is  provided  here. 
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A more  complete  discussion  of  the  numerical  procedure  used  to 

obtain  da!,  and  dP  under  conditions  of  plastic  yielding,  is 

1J  A6 
found  elsewhere 

According  to  the  method  of  the  plastic  potential,  the 
incremental  stress-strain  relation  during  shear  failure  — a 
"flow  rule"  — can  be  completely  determined  if  the  locus  of 
stress  states  in  which  shear  failure  occurs  — the  "yield  sur- 
face" — is  known  for  a given  material.  Here  the  yield  surface 
is  described  as  a combination  of  a generalized  Mohr-Coulomb 
surface,  and  a von  Mises'  surface.  The  yield  criterion  can  be 
written  in  form: 

/Sj  * Y(P)  (A15) 

where  J^,  the  second  invariant  of  the  deviatoric  stress  ten- 
sor, is  defined  by  the  equation 

J2  = Kj  aij  <A16> 

The  "yield  function",  Y (P)  , is  defined  by  the  relation 


Y(P) 


P<P 


P^P 


(l-e/e)(l-(e/e)H) 
m v 


(A17) 


where  Y^M  and  YMC  are  yield  functions  of  von  Mises  and  gen- 
eralized Mohr-Coulomb-type,  respectively,  Pq  is  the  mean  stress 
at  which  the  transition  from  Mohr-Coulomb  to  von  Mises'  yield 

A 

functions  occurs,  and  the  factors  (1-e/e  ) and  [l-(e/e  ) ] ac- 
count  for  thermal  weakening  of  the  material.  Figure  A2  contains 
a schematic  of  the  yield  function,  Y(P),  used  in  the  present  pro- 
gram. 
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Figure  A2:  Typical  Mean  Stress  (?)  Dependent  Yield  Function  (Y(P)).  Shear 
Failure  occurs  when  the  square  root  of  the  second  invariant  of  the  deviatoric 
stress  tensor,  \/jX,  is  equal  to  Y(P). 


A2.3  Calculation  of  Stresses  during  Shear  Failure:  The 
Flow  Rule 

Application  of  stress  to  an  element  of  material  ("loading") 
causes  deformation  of  the  element.  Removal  of  the  stress  ("un- 
loading") generally  results  in  only  partial  return  of  the  ele- 
ment to  its  undisturbed  configuration.  The  part  of  the  con- 
figuration change  that  is  reversed  or  "recovered"  on  unloading 
is  termed  the  "elastic"  part  of  the  deformation  caused  by  load- 
ing, and  is  denoted  De<  The  net  overall  change  in  configuration 
after  loading  and  unloading  is  termed  the  "inelastic"  (or  "plastic") 
part  of  the  deformation  caused  by  loading,  and  is  denoted  D^. 

The  configuration  change  produced  on  loading  can  thus  be  viewed 

as  the  result  of  the  inelastic  deformation  D on  the  undisturbed 

P 

element  followed  by  the  elastic  deformation  D . 

e 

The  strain  tensor,  e„,  provides  a quantitative  descrip- 
tion of  strains  associated  with  the  elastic  and  inelastic  parts 

of  that  deformation  denoted  by  and  e?}.  The  total  strain 

el  1Jpl  1J 

e..  is  simply  the  sum  of  e . . and  eV.:  in  incremental  form,  we 
ij  r J ij 

write 


de . . = de?l  + dc 
ij  iJ  iJ 


(A18) 


The  plastic  strain  increment,  de?l  is  related  to  the 

ij 

stress  increment  through  the  flow  rule 


depl  = *£ 
ij  9a 


dx 


ij 


F.  .dX 
IJ 


(A19) 


where  F is  the  plastic  potential.  Here  an  "associated"  flow 
rule  has  been  adopted;  the  plastic  potential  F is  defined  by 
the  yield  surface: 
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F = J ' 
J2 


Y (P)  = 0 


(A20) 


By  combining  Eqs . (A13) , (A18) , (A19) , and  (A20) , the  follow- 
ing hydrostatic  and  deviatoric  incremental  stresi  strain  re- 
lations for  plastic  flow  are  obtained: 


del.  = 


da!  . 

-sr1  + F! .dX 
2G  ij 


dA  = - + cp  dX  = dAe^  4-  cpdX 


(A21) 

(A22) 


,dY 


where  ^F.^-ZY^  and  F'.  = F..  - 6..  cp/ 3 = a!.. 


ij 


ij 


ij 


By  contracting  Eq.  (A21)  with  aj^  to  obtain 
dX  = [ o[.  de!  . + c'jda[j/2G]/2J^ 

and  recognizing  that  a! .da' . = dJi  = 2Y-—  and  dF  = -kdAel, 

lj  lj  z dp  ’ 

Eqs.  (A21)  and  (A22)  can  be  written  in  the  form 


da!  . = 
ij 


■2G[dEij  - 


+ »i3Y(p>§  • ^ 


ei 


(A23) 


(A24) 


, dA  + 4^(a!  de ! . ) /Y 
dAel  

[1  + -(^X)  ] 

L GMP;  J 


(A25) 


where  k is  the  bulk  modulus,  G is  the  shear  modulus  and  Y(P) 
is  defined  by  Eq.  (A17) . 

For  uniaxial  and  spherically-symmetric  motion,  the  stress 


K< 


i 
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axes  have  known  fixed  directions,  and  two  of  the  principal 
stresses  are  equal.  In  addition,  the  principal  stress  devia- 
tors and  strain  deviator  increments  satisfy  relations  of  the 
form  The  general  incremental  stress-strain 

equation  (A25)  then  reduces  to  the  simpler  form: 


where 


dA 


el  _ + sgn<°x)v'Tl?  K 

1 + (— ) (— ) 


Sgn(°x)  = 'x'K 


(A26) 


(A27) 


Y = 


(A28) 


and  x designates  the  direction  of  motion.  To  complete  the 
constitutive  model,  the  sequence  of  calculations  must  be  de- 
fined more  precisely;  a few  computational  details  that  be- 
come important  in  the  event  of  shear  failure,  also  need  to 
be  stated. 

A3.  The  Steps  Taken  to  Compute  Stress  Increments  from  Strain 

Increments 

The  first  step  in  the  series  of  calculations  whereby  any 
stress  increment  is  computed,  is  to  determine  whether  the  de- 
viatoric  strain  increment  is  elastic  or  not.  To  answer  that 
question,  the  deviatoric  strain  increment  is  assumed  at  the 
outset  to  be  elastic.  On  that  assumption,  Eq.  (A13)  and  the 
mean-stress  equations  of  Section  A1  are  used  to  calculate  a 
virtual  stress  increment  and  virtual  stress.  If  the  virtual 
updated  stress  satisfies  the  inequality  (A15),  then  the  virtual 
increment  and  state  are  taken  as  actual.  Otherwise,  the  devia- 


A14 


toric  strain  increment  is  not  elastic,  and  Eqs.  (A13)  through 
(A18)  and  (A24)  and  (A25)  are  used  to  compute  the  actual  stress 
state. 

If  the  deviatoric  strain  increment  is  not  elastic,  then 
the  elastic  part  of  the  dilatation  increment  is  computed  ac- 
cording to  Eq.  (A25).  The  bulk  modulus  appropriate  to  the 
virtual  state  of  updated  stress  is  used  for  that  purpose  (ex- 
perience has  shown  that  this  easily- improved  approximation  is 
adequate).  The  excess  compression,  p,,  that  enters  the  hydro- 
static stress-strain  relations  is  then  updated  using  only  the 
elastic  part  of  the  dilatation  increment. 

A4.  Equation  of  State  for  Wet  Sandstone 

The  equation  of  state  for  wet  sandstone  was  developed 

from  basic  material  properties  data  from  several  sandstones. 

A7 

Wavespeed  and  density  data,  provided  to  us  by  the  USGS  , are 
for  sandstones  from  the  region  of  interest  in  the  Soviet  Union. 
Other  data  were  obtained  from  U.S.  sandstone  and  silica  and 
were  used  to  develop  the  general  equation  of  state.  With  the 
parameter  values  listed  in  Table  Al,  Eqs.  (A3),  (A7),  (A2)  and 
(A14)  define  ev,  p,',  and  G respectively.  Mean  stress  for 
the  solid  phase  is  described  by 


V' 

p'^O 

F(H\  »*,  Ko,  KJ 

o 

A 

The  definition  of  the  shear  modulus  is  completed  with 

G(u  ) * MINCG  ,G  + P_(n'  _)] 

max  m o dr  s max 

For  this  saturated  material  the  yield  criterion  is  independent 


(A29) 


(A30) 
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TABLE  Al 


MATERIAL  PARAMETERS  FOR 
SOVIET  GRANITE,  AND 


WET  SANDSTONE 
NTS  GRANITE 


Variable  (Units) 


Wet 

Sandstone 

Soviet 

Granite 

NTS 

Granite 

2.400 

2.700 

2.661 

0.416666666 

0.370370370 

0.375798572 

3.050 

2.980 

2.000 

69.120 

522.500 

256.000 

550.000 

665.000 

940.000 

0.050 

0.015 

0.0279 

-25.000 

-40.000 

-40.000 

0.000 

0.000 

0.000 

1.000 

1.150 

3.500 

3.500 

3.500 

3.500 

2.115 

5.169 

0.000 

5.000 

5.360 

3.500 

0.000 

0.000 

0.000 

5.000 

5.000 

5.000 

0.100 

0.100 

0.100 

0.400 

0.400 

0.400 

0.900 

1.300 

1.300 

10.000 

16.000 

16.000 

51.840 

194.000 

I 

- 

120.000 

339.500 

- 

11.7885 

- 

- 

- 

0.005 

- 

250.000 

577.350269 

3000.000 

of  the  mean  stress  P.  The  Mohr-Coulomb  surface  is  removed  so 
that  for  all  P 


Y(P)  = Yvm  (A31) 

Plots  of  material  properties  for  the  wet  sandstone  model  are 
given  in  Figures  A3  through  A7. 

A5.  Equation  of  State  for  Soviet  Granite 

As  with  wet  sandstone  the  material  model  for  Soviet 
granite  is  based  on  properties  provided  by  USGS  for  several 
granites.  Data  available  for  Soviet  granite  included  longi- 
tudinal wavespeed  as  a function  of  pressure,  a more  detailed 
picture  than  that  available  for  typical  Soviet  sandstones.  The 
overall  model  is  based  on  the  thermodynamic  properties  of 
several  granitic  rocks.  The  equations  used  to  define  this  wet 
granite  (non-hysteretic)  are  the  same  as  those  used  for  wet 
sandstone  except  that  the  shear  modulus  is  given  by 


% 

■ 


* 


G (p,  ) = Max[0,  G -(G  -G  )exp(-p  /p,  )]  (A32) 

’max'  * m v m o'  max  g 

and  parameter  values  are  those  for  Soviet  granite  given  in 
Table  Al.  Characteristic  properties  are  plotted  in  Figs.  A8 
through  A12. 

A6.  Equation  of  State  for  Dry  Sandstone 

The  model  for  dry  sandstone  is  the  most  complicated  model 
derived  for  this  program.  The  data  base  was  the  same  as  that 
provided  ATI  by  the  USGS  for  wet  sandstone.  Neither  the  pres- 
sures at  which  cementation  breakdown  and  pore  collapse  occur 
nor  data  on  the  strength  of  this  dry  material  were  provided.  On 
the  basis  of  discussions  with  personnel  from  USGS  and  Terra  Tek, 
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Figure  A4.  Properties  of  Wet  Sandstone:  Shock  Speed  S and  Shear  Wave 
Speed  Cr  vs.  Stress  for  Uniaxial  Strain  Conditions. 
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Figure  A12. 


Properties  of  Soviet  Granite:  Hugonlot  and  Release  Adlabat 
Stress  (o  ) vs.  Compression  (T])-High  Stress  (q  > lOOkbar) 


data  for  the  sandstone  from  the  sites  of  the  Wagon  Wheel  or  Gas 
Buggy  Events  were  considered  appropriate  for  defining  these 
features  of  the  dry  Soviet  Sandstone  model. 


The  equations  for  the  vapor  portion  of  the  model  are 

the  same  as  those  in  the  wet  sandstone  equation  of  state,  as 

are  the  equations  for  the  energy  dependence  of  G and  Y.  The 

equations  for  the  loading  and  unloading  hvdrostats,  P , are 

s 

given  below.  These  differ  significantly  from  those  for  the 
other  media  treated  here  due  to  the  complex  changes  occuring 
during  cementation  breakdown.  The  loading  hydrostat  ) 

is  given  by 


and 


[ Ki  p' 

> P ' 

£ 

^e 

Ps  = j 

Sp(p, i) 

; ^e 

< 

u'<p  ; 1 < i < N 
m 

; p' 

2 

P 

m 

p 

max 


MIN(p\  nm) 


where  Sp(p,',i)  are  cubic  splines  describing  the  P-^  curve  over 
N ( = 14)  regions  of  the  excess  compression  interval  p,  < y, ' < 
The  unloading  hydrostat  ) is  given  by 


P 
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V 

Pu(“">lVax) 

F(^",u*,KQ.Km) 


u ^ \x 

max  e 

U < til  1 

e max 

p. ' s p 

max  m 
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P 
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The  function  pu(p" >P^ax)  represents  the  complex  inter- 
mediate range  of  the  unload-reload  hydrostat.  Its  functional 
form  is  that  of  Eq.  (A9)  , but  the  following  series  of  steps  is 
necessary  to  define  its  arguments  and  guarantee  P,^  compati- 
bility at  the  load-unload  transition. 


(A33) 

(A34) 


(A35) 


Pmax  “ SP^max* 


^ = S (p'  , i) 

z p v ^max  * 


(A36) 


(A3  7) 


where  Sp(p' ,i)  is  the  cubic  spline  representing  the  appropriate 

segment  of  the  loading  hydrostat  and  S Is  the  related 

cubic  spline  representation  of  the  permanent  excess  compression 

p that  will  occur  on  unloading  from  p'  . From  p the  variable 
z max  z 


i"  is  defined. 


• t « / I \ 

IX  “ LL  - U (u,  ) 

max  max  z max 


The  value  for  Kl  is  then  found  from  the  equation 


P = F(p"  , p*,  K',  K ) 
max  max’  * o*  m 


The  pressure  is  then  given  by 


(A38) 


(A39) 


Pu  - F(,",  u*,  k;,  Km) 


(A40) 


The  shear  modulus  for  dry  sandstone  is  given  by 


G(p  ) = Min(Qf£  ,G  ) 

v>max  op  , , * m 

unload 


(A41) 


Finally,  the  yield  surface  is  defined  by  the  relation 


Y1  + Y2  P 

Y (P)  = |y3  + (Yvm-Y3)(1-(1-Z)4) 


P s P. 


P < P < P 
MC  \ 


P * P, 


(A42) 


where  Z = — - . This  form  produces  a smooth  transition  from 

VM  MC 

the  Mohr-Coulomb  form  to  the  von  Mises  limit  over  the  interval 


t 
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PMC  < P < PVM  as  shown  in  Figure  A2. 


Values  of  the  parameters  appearing  in  these  expressions  are 
given  in  Table  A2.  Coefficients  for  the  cubic  splines  defining  Sp 
and  S are  given  in  Tables  A3  and  A4,  respectively.  Plots  of  char- 
acteristic properties  are  shown  in  Figures  A13  through  A19. 

The  model  for  weak  dry  sandstone  is  the  same  as  that  for 
dry  sandstone  except  that  the  yield  surface  is  that  of  the  wet 
sandstone  model  (Section  A4) . Plots  of  the  material  properties 
for  this  hybrid  model  are  given  in  Figures  A20  through  A23. 

A7.  Equations  of  State  for  NTS  Granite 

A1 

The  model  for  NTS  granite  was  developed  in  the  late  1960s 
during  a period  when  the  Piledriver-Hardhat  dilemma  was  being 
studied.  The  model  for  the  high  energy  behavior  of  the  material 
was  substantially  more  complicated  than  the  current  model  although 
the  basic  equations  are  the  same  as  those  used  here  for  Wet  Sand- 
stone and  Soviet  granite.  The  equations  defining  the  shear  be- 
havior for  NTS  granite  differ  from  those  for  the  other  media. 

The  shear  modulus  is  given  by: 


(A43) 


where  the  constant  0.6  implies  a constant  Poisson's  ratio 
(v*0.25).  The  yield  surface  is  defined  by 


where 


Y (P)  = MIN  (Y  , Y +Y,  P) 
v ' ' max’  o 1 ' 


Y = Y,_,(A  - Be) 
max  VM 


(A44) 


(A45) 


Yq  - 100  bars,  Y^  - 0.500,  A=1.67,  B=0. 67x10  gm/erg  and 


Yvm  - 3000  bars. 
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TABLE  A2 

MATERIAL  PARAMETERS  FOR  DRY  SANDSTONE 


SPLINE  COEFFICIENTS  FOR  DRY  SANDSTONE  LOAD  HYDROSTAT 
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Figure  A16.  Properties  of  Dry  Sandstone:  Permanent  Excess  Compression, 
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Figure  A17.  Properties  of  Dry  Sandstone:  Bulk  Modulus  K,  Shear  Modulus  G 
Poisson's  Ratio  v,  P-Wave  Speed  C and  S-Wave  Speed  C vs. 
Mean  Stress  P for  Hydrostatic  Loading. 
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The  values  of  parameters  appearing  in  the  equations  for 
the  hydrostat  are  listed  in  Table  A1  with  those  for  wet  sand- 
stone and  Soviet  granite.  Plots  of  characteristic  properties 
are  given  in  Figures  A24  through  A28. 
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Figure  A25.  Properties  of  NTS  Granite:  Shock  Speed  S and  Shear  Wave 


Speed  Cg  vs.  Stress  for  Uniaxial  Strain  Conditions. 


Figure  A26.  Properties  of  NTS  Granite:  Bulk  Modulus  K,  Shear  Modulus  G 
Poisson's  Ratio  v,  P-Wave  Speed  C.,,  and  S-Wave  Speed  Cg  vs. 
Mean  Stress  P for  Hydrostatic  Loading. 
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Figure  A28.  Properties  of  NTS  Granite;  Hugoniot  and  Release  Adiabat  Stress 
(a  ) vs.  Compression  (n)  -High  Stress  (ax  > lOOkbar)  Region. 
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APPENDIX  B 
BODY-WAVE  MAGNITUDES 

Bl.  Body-Wave  Magnitudes  at  Teleseismic  Distances:  A Simple 

Model 

Approximate  values  of  body-wave  magnitude,  m^,  are  computed 
using  a procedure  which  is  exact  to  the  extent  that  explosions  in 
the  earth  generate  spherically  symmetric  fields  in  linear  media. 

Ground  motions  from  buried  explosions,  however,  differ 
from  the  idealized  spherically  symmetric  fields  in  linear  media 
in  two  obvious  ways.  First,  near  the  explosive  source  the  motion 
is  so  violent  that  strong  nonlinear  behavior  occurs.  Second,  for 
problems  of  practical  interest  the  ground  surface  disturbs  the 
approximate  spherical  symmetry  producing  more  complex  2-dimensional 
axisymmetric  fields  of  motion.  These  deviations  from  the  idealized 
behavior  used  in  computing  body-wave  magnitudes  can,  to  some  ex- 
tent, be  accounted  for. 

Nonlinear  inelastic  behavior  will  occur  out  to  some  medium 
and  burst-specific  slant  range  (spherical  radius)  from  the  shot 
point.  Beyond  that  range  only  elastic  deformation  will  occur. 

By  monitoring  computed  ground  response  in  an  evolving  field  of 
motion,  the  maximum  range  for  inelastic  deformation  can  be  deter- 
mined. For  points  at  larger  ranges,  pulses  of  particle  velocity, 
stress,  and  other  response  variables  then  truly  represent  linear 
(elastic)  response. 

For  one-dimensional  spherically  symmetric  motion  no  sym- 
metry approximation  is  involved.  Computed  velocity  pulses  ob- 
tained at  several  ranges  as  indicated  in  Figure  Bl(a)  can  be 
used  directly  to  describe  the  field  of  motion.  For  two-dimen- 
sional axisymmetric  motion,  vertical  and  horizontal  components 


Bl 


a)  Spherically  Symmetric  b)  Axisymmetric  2-D  Fields 
One-D  Fields 


Figure  B-l.  Velocity  Pulses,  UR,  for  Computing  Body-Wave 
Magnitudes . 


of  particle  velocity  and  other  data  were  obtained  for  points  at 
1°  intervals  along  the  90°  arc  in  the  quarter  plane  (flow  plane) 
shown  in  Figure  Bl(b).  For  each  station  considered,  the  vector 
component  along  the  radius  OP  from  the  shot  point  was  taken  as 
an  equivalent  spherical  particle  velocity  pulse.  This  pulse 
varies  from  point  to  point  reflecting  differences  in  arrival 
times  and  amplitudes  of  contaminating  signals  from  the  ground 
surface.  These  equivalent  spherical  pulses  from  stations  out- 
side the  inelastic  regime  were  the  input  data  for  computing 
body-wave  magnitude. 

Once  velocity  pulse  data  is  available,  our  computation 
of  m^  consists  of  four  steps  (Section  B2) : 

(1)  Treat  the  radial  velocity  pulse  computed  for  a point 

(slant-range  r^)  where  deformation  is  purely  elas- 
tic as  an  exact  representation  of  a spherical  P-wave. 
With  the  pulse  at  r^  as  boundary  condition,  use  the 
equations  of  continuum  motion  to  calculate  the  pulse 
that  would  be  observed  in  an  infinite  elastic  medium 
at  a much  greater  (teleseismic)  range  of  interest 
(range  r2). 

(2)  Compute  the  Fourier  spectrum  of  the  pulse  at  r2» 

(3)  Damp  each  Fourier  component  by  the  factor  exp  (-u>t/2Q)  , 
where  w is  frequency,  t is  the  time  of  travel  of  a P- 
wave  from  r^  to  r2,  and  Q is  a constant  (1000  in  most 
of  the  calculations  reported. 

(4)  Superpose  all  damped  Fourier  components  to  produce  an 
attenuated  pulse  at  r2> 

The  pulse  from  step  (4)  is  used  as  a forcing  function  to  drive 
a mathematical  seismometer  at  r2;  m^  is  computed  from  the  maxi- 
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mum  peak-to-peak  displacement  of  the  resulting  seismometer  re- 
sponse (Section  B5). 

In  the  bare-bones  procedure  outlined  above,  two  major 
complications  have  been  bypassed:  (a)  ref lection- transmission 
from  an  interface  between  the  local  near-field  medium  and  the 
rest  of  the  earth,  and  (b)  phase  and  amplitude  corrections  made 
in  step  (4)  so  that  motion  begins  at  r ^ at  the  time  when  a P- 
wave  first  arrives  at  that  range,  and  not  before. 

The  need  for  item  (a)  rises  from  treating  the  earth  as  a 
standardized  vehicle  for  seismic  wave  propagation  to  provide  a 
reasonable  basis  for  comparing  m^-values  for  bursts  in  dif- 
ferent local  media.  To  reconcile  the  differences  in  local  media 
with  the  standard  earth  a spherical  interface  was  created  at  dis- 
tance rp  from  the  shot  point,  where  r^<rp«r2.  The  transmission 
coefficient  for  wave  propagation  across  the  interface  (Section 
B4)  was  obtained  by  solving  the  equations  that  govern  spheri- 
cally symmetric  continuum  motion.  (r  =5  km  for  the  calculations 

r 

reported,  while  at  greater  ranges  the  earth  was  a homogeneous 
linear  elastic  medium  with  P-  and  S-wave  speeds  of  8.0  km/sec 
and  4.6  km/sec,  respectively. 

Relative  to  item  (b) , it  is  well  known  that  if  step  (4)  is 

followed  literally,  then  disturbances  will  reach  teleseismic 

ranges  instantaneously.  Incomplete  experimental  definition  of 

the  earth's  Q-factor  is  responsible  for  that  causality  violation 

[the  same  problem  is  met  when  electromagnetic  waves  travel  through 

a medium  whose  dielectric  constant  is  known  only  on  a limited 

Bl 

frequency-range].  In  our  m^  calculation,  a standard  method  is 
used  to  recover  causality,  while  sti’l  reproducing  measured  dis- 
persive properties  (Section  B3). 


B4 


jmmMA 


While  important,  refinements  (a)  and  (b)  do  not  lie  at 
the  heart  of  our  calculation  of  m^;  the  basic  concept  of  the 
procedure  is  the  propagation  of  spherical  waves  through  a dis- 
sipative medium.  A more  detailed  discussion  of  the  method, 
including  corrections  (a)  and  (b) , is  given  in  the  following 
sections . 

B2.  Homogeneous  Isotropic  Medium:  No  Causality  Correction 

Let  v(r,t)  be  the  radial  particle  velocity  at  time  t and 
at  radius  r from  a center  of  spherical  symmetry,  in  a homoge- 
neous, isotropic,  linear  elastic  medium  whose  P-wave  speed  is 
c.  Then  v(r1}t)  denotes  particle  velocity  as  a function  of 
time  at  radius  ra . Given  v(r!,t),  an  exact  explicit  expres- 

B2 

sion  can  be  written  for  v(r,t);  simple  quadrature  suffices 
to  generate  the  particle  velocity  at  any  point  and  time,  and 
in  particular  at  radius  ra=r1+Ar»r1,  as  a function  of  time. 

Decomposing  a velocity  or  displacement  field  into  Fourier 
components,  and  multiplying  the  component  of  frequency  uu  by 
exp(-  (Ar/2Qc) , is  perhaps  the  simplest  (and  most  widely  used) 
method  of  accounting  for  the  fact  that  waves  are  damped  as 
they  travel  through  the  earth.  In  the  present  case,  a damped 
pulse  v*(ra,t)  is  obtained  by  adding  the  resulting  harmonics; 
v*(r3,t)  more  nearly  represents  the  input  to  a seismometer  at 
ra  than  v(ra ,t) . 

In  presenting  a detailed,  self-contained  account  of  the 
calculation  of  v*(r2,t)  here,  we  first  write  the  equations 
that  govern  v(r,t)  and  its  Fourier  transform,  v(r,u>): 
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r dr 
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where 


v(r  ,u>)  = 


+® 

y v(r. 


t)e1(,tdt 


Evidently,  in  the  pulse  that  arrives  at  radius  ra , v(ra,w)  is 
the  Fourier  component  of  frequency  w . Hence,  taking  account 
of  damping  by  the  earth,  the  pulse  felt  by  a seismometer  at 
ra  is  given  by 


v*(r. 


A,  s -|iulAr/2cQ  -iwt, 
v(r3  ,u>)e  x e d^ 


To  evaluate  the  integral  of  Eq.(B4),  v(ra  ,w)  is  needed. 
Eq.(B2)  is  used  to  obtain  v(ra,ui).  Since  all  sources  of  radia- 
tion lie  within  a sphere  of  radius  r1«ra  , only  an  outgoing 
wave  is  permitted  at  each  frequency;  Eq.(B2)  then  implies  that 

0(r,»)  - A(.)  i2(l  +i£)e1“r/c 

The  amplitude  A (w)  is  found  using  the  known  particle-velocity 
pulse  at  rx . In  particular,  it  follows  from  Eqs.(B3)  and  (B5) 
that 
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Eqs . (B4)- (B6)  yield  directly  the  following  equation  for  v*(ra,t): 


+ oo 


+ ® 


v*(r»»c)  = 


v(r1,t')e1U,t  dt', 


|“>  |Ar/ 2cQ 


B((u)e  ^U)tdu)  (B7) 


where 

B(»)  - ei*Ar/c  (B8) 

r3  (l+ic/wr!) 

Eq.(B7)  is  exact  for  a spherical  field  in  a dissipative, 
but  otherwise  isotropic,  linear  elastic  medium  - provided  that 
dissipation  can  be  represented  for  any  frequency  u>  by  the 
factor  exp(- |u>  |Ar/2cQ)  . 


Two  important  modifications  to  Eq.(B7)  are  made  below. 
Before  turning  to  those  modifications  in  detail,  we  observe 
that  by  associating  the  exponential  damping  factor  with  the 
inner  integral  of  Eq.(B7),  the  equation  can  be  interpreted  as 
describing  linear  elastic  propagation  to  radius  ra , of  a damped 
pulse  at  rj . Since  waves  actually  damp  as  they  travel,  that 
interpretation  is  conceptually  unattractive.  However,  the 
mathematical  equivalence  of  the  two  viewpoints  is  important: 
it  makes  absolutely  no  difference  whether  (a)  the  near-field 
pulse,  modified  to  account  for  damping,  propagates  to  the  far 
field  in  perfectly  elastic  fashion,  or  (b)  the  near-field 


pulse  propagates  to  the  far-field  in  perfectly  elastic  fashion 
and  is  then  damped.  The  same  is  true  in  modifying  Eq.(B7)  both 
to  take  account  of  near-field  reflection,  and  to  preserve 
causality  in  the  far  field.  Each  of  those  two  modifications 
gives  rise  to  an  additional  factor  in  the  integrand  of  Eq.(B7) 
and  that  factor  can  be  considered  as  altering  directly  either 
the  near-  or  far-field  pulse.  Physically,  of  course,  reflec- 
tion from  a near-field  interface  modifies  the  near-field  pulse, 
and  a causality  correction  is  needed  not  in  the  near  field  but 
in  the  far  field. 

B3.  Causal  Wave  Propagation  with  Damping 

With  the  wavespeed  c constant,  causal  propagation  to  an 
arbitrary  radial  position  r>rT  is  obtained,  together  with 
damping,  from  the  formula 


v(r,u>)  = v(rT  ,w)B(r  ,r!  ,w;c)H(r,r1  ,ui;c) 


(B9) 


where 


B(r ,rx  , uj ; c ) 


/ 1 + ic/wr  \ iw(r-rl)/c 
r \ 1 + ic/ mr, j e 


(BIO) 


H(r ,ri ,w;c)  = e" * ” L*i 


(Bll) 


The  variables  * and  ^ (Eq.(Bll))  are  defined  as  follows 


B1 


♦ 


lll£i 1 

2Qc 


min(  | id  | ,wH) 


(B12) 


B8 


where 


g H lUJH/tul  (B14) 

Up  to  now,  wH  has  been  set  at  5000  rad/sec,  and  Q is  presently 
taken  as  1000;  also  500m<r1 <2500m.  [Note:  the  formulas  above 
may  be  wrong  for  that  case  has  yet  to  arise  in  numerical 

practice . ] 

B4.  Causal  Propagation  to  Teleseismic  Distances,  Including  Re- 
flection from  a Near-Field  Interface:  Present  ATI  Equations 

We  view  the  earth  as  one  and  the  same  homogeneous  medium 

over  almost  the  entire  length  of  any  teleseismic  propagation- 
path.  At  the  same  time,  the  mechanical  properties  of  geologic 
media -exhibit  wide  local  variations.  To  account  for  both 
local  variability  and  teleseismic  constancy  of  the  earth, 
the  effects  of  an  interface  between  the  two  media  - local  and 
global  - are  included  in  our  equations.  The  interface  appears 
at  radius  rp  from  a given  shot  point  (in  most  calculations  to 
date,  rp=5000m) . Relative  to  passage  of  signals  across  the 
interface,  both  the  local  and  global  media  are  considered 
homogeneous,  isotropic  linear  elastic  materials. 

Subscripts  "1"  and  "a"  are  used  below  to  refer,  respec- 
tively, to  properties  of  the  local  medium  (r<r_)  and  surround- 

r 


ing  earth  (r>r  ) . 

r 

In  computing  particle  velocities  at  teleseismic  distances, 
the  pulse  incident  upon  the  interface  is  first  obtained  from 
Eqs . (B9)- (B14) . For  that  purpose  r and  C are  replaced  by  rf  and 
Cj  in  the  equations  cited.  The  resulting  component  of  interface 
velocity,  vCr^jU)),  then  satisfies  the  relation 

v(rF,ui)  = v(rx  ,uj)B(rF,ra  , w;cT )H(rp ,r*  ,u>;Cl  )T(u>)  (B15) 

where  the  transmission  coefficient  T(w)  is  given  by  the  expres- 
sion 
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In  Eq.  (B16),  ax  and  ag  are  shear-wave  speeds.  Also,  at  present, 
Pa-2. 7 gm/cm3,  c*  = 8 km/sec  and  ae=S/^3  km/sec;  pT  , cx  and  ax 
vary  with  the  local  medium. 

To  complete  the  calculation  of  velocity  at  a teleseismic 
range  ra , a function  G(t)  is  defined  as  follows: 


G(t) 


+“ 


(B17) 
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where  *'  and  replace  t and  * respectively,  in  Eqs . (B12)- (B14) 
and  are  obtained  from  those  equations  by  substituting  r for  rj , 
ra  r>  a°d  c8  for  c.  The  function  G(t)  and  the  particle 
velocity  v(ra,t)  at  r8 , are  related  by  the  equation 


rF 

<L  G/t 

- r°"rF)  + £tG/t 

rB'rF\l 

ra 

dtGV 

c.  / + 7TG1‘ 

ca  / 

Eq.(B18),  in  the  following  discrete  approximation,  is  used  to 
compute  v(ra , t) : 


v(ra , t+At 


where 


e = (ca/rF)^t  ( 

@2  = (ce  /r2  )At  ( 

All  Fourier  Transform  integrals  are  evaluated  by  the 
Fast  Fourier  Transform  method. 

B5.  Seismometer  Response 

The  velocity  history  at  r ^ which  we  have  obtained  by  the 
procedure  described  above  is  used  to  drive  a mathematical  seis- 
mometer. The  response  of  the  seismometer  is  then  used  to  obtain 
the  body-wave  magnitude.  Representations  of  two  seismometer 

models  have  been  used  during  this  program.  The  first  model  (pro- 

R 1 

vided  ATI  by  DARPA/NMRO)  is  that  of  a WWSSN  seismometer  . The 
Response  Function  of  the  frequency  dependent  model  is  given  by 
the  complex  analytical  expression: 


(B19) 

(B20) 
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S = 18 
r = 170 
s = 18 
and  i =* 


The  second  seismometer  model  is  for  an  AFTAC  instrument.  This 

model  is  the  same  as  the  one  used  by  System,  Science  and  Soft- 

84 

ware  in  their  calculations  of  magnitudes  . In  this  model  the 
response  of  the  instrument  is  also  a complex  function  of  the 
angular  frequency,  where 
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R.)  - I(»/2n)e1,p("/2n) 

The  amplitude  and  phase  of  this  frequency  dependent  model  are 
obtained  by  interpolating  in  a table  of  tp(v)  and  log1QI(v) 
where  ®(v)  is  the  phase  angle,  I(v)  is  the  amplitude  of  the 
seismometer  response  and  v is  frequency.  Table  B1  gives  the 
amplitude  and  phase  at  discrete  frequencies. 

To  obtain  the  response  of  the  seismometer  to  the  input 

A 

velocity  at  r?,  the  velocity  is  Fourier  transformed  to  VCr^jUj). 
The  spectrum  of  the  displacement  at  r ^ is  calculated  from  the 
velocity  by  dividing  the  velocity  spectrum  by  i<u.  The  displace- 
ment spectrum  is  then  multiplied  by  the  response  function  of  the 
instrument  to  give  the  spectrum  of  the  seismometer  displacement, 

6(ti)) 

6(u>)  = n>)[V(r2,u>)/iuj] 

The  spectrum  is  then  transformed  into  the  time  domain  to  give 
the  seismometer  deflection  6 (t ) . 

To  determine  the  magnitude  of  an  event  from  the  seismom- 
eter record  the  following  equation  is  used: 

= logl0(A/(TI(T))  + B(A) 

where  A is  the  maximum  peak  to  trough  amplitude,  T is  twice  the 
time  between  the  peak  and  the  trough,  I(T)  is  an  instrument 
correction  factor  which  is  a function  of  the  period,  T,  and  B(A) 
is  a correction  factor  for  the  distance  from  the  location  of  the 
event  to  the  seismometer  station.  For  all  magnitudes  reported 
here,  the  distance  is  5000  km  (A =45°)  and  the  correction  B(A) 
is  3 . 2885 . 
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TABLE  B1 


AMPLITUDE  AND  PHASE  OF  AFTAC 
SEISMOMETER  RESPONSE  FUNCTION 


Index 

v(hz) 

I 

cp(rad) 

1 

.0100 

.0000001 

6.10863 

2 

.0250 

.0000045 

5.58503 

3 

.0500 

.000075 

5.20106 

4 

.1000 

.001 

4.77869 

5 

.2000 

.010 

4.40519 

6 

.4000 

.080 

3.83971 

7 

.5000 

.145 

3.65820 

8 

.6000 

.243 

3.42083 

9 

.7000 

.375 

3.19394 

10 

.8000 

.548 

3.00209 

11 

.9000 

.750 

2.80997 

12 

1.0000 

1.000 

2.55865 

13 

1.2500 

1.650 

1.89367 

14 

1,5000 

2.480 

1.36484 

15 

1.7500 

2.890 

0.76794 

16 

2.0000 

3.030 

0.31939 

17 

2.2500 

2.980 

-0.05236 

18 

2.5000 

2.860 

-0.39968 

19 

2.7500 

2.720 

-0.68068 

20 

3.0000 

2.580 

-0.93200 

21 

3.5000 

2.290 

-1.39626 

22 

4.0000 

2.040 

-1.75929 

23 

4.5000 

1.800 

-2.13802 

24 

5.0000 

1.570 

-2.43123 

25 

6.0000 

1.190 

-3.00719 

26 

7.0000 

.880 

-3.51508 

27 

8.0000 

.643 

-3.96712 

28 

10.0000 

.336 

-4.73157 

Figures  B2  and  B3  are  plots  of  I versus  v and  the  pro- 
duct TI  versus  T for  the  AFTAC  seismometer  while  Figures  B4 
and  B5  give  similar  plots  for  the  WWSSN  seismometer. 
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Figure  B2 


Amplitude  of  the  Instrument  Response  as  a Function  of  the 
Frequency,  for  the  AFTAC  Seismometer. 


FREQUENCY  (HZ) 


Figure  B4 . Amplitude  of  the  Instrument  Response  as  a Function  of  the 
Frequency,  v,  for  the  WWSSN  Seismometer, 
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APPENDIX  C 
MODEL  DESCRIPTION 

Cl.  Axisymmetric  Source  Calculations  (2D)  for  Buried 

Explosions 

The  matrix  of  calculations  adopted  for  this  program 
shortly  after  it  began  consisted  of  twelve  bursts  in  four  dif- 
ferent media.  As  the  program  evolved,  the  problem-set  was 
expanded  to  25  bursts  in  five  media,  as  specified  in  Table  Cl; 
the  geologic  materials  and  depths  of  burial  for  all  problems 
are  listed  therein.  The  general  method  of  calculation  (Sec- 
tion C2.)  employed  for  the  original  12  problems,  difficulties 
discovered  in  using  the  output  from  those  calculations,  and  a 
revised  method  of  calculation  which  avoids  those  difficulties 
are  described  in  the  following  subsections. 

C2.  Original  Computational  Technique 

Computing  explosively  driven  ground  motion  is  not  just 
a simple  matter  of  feeding  into  a computer  initial  conditions, 
material  properties,  and  a finite  difference  mesh,  and  waiting 
for  the  answers.  To  achieve  desired  levels  of  accuracy  in 
describing  ground  response  while  meeting  problem  size,  running 
time,  and  cost  limitations  typically  requires  that  the  calcula- 
tion be  performed  in  several  stages.  During  each  stage  the 
mesh  can  be  chosen  to  treat  specific  features  critical  to  that 
stage,  e.g.,  rock  vaporization  and  cavity  expansion  during  the 
initial  stage,  or  ejecta  production  during  an  intermediate 
stage.  Details  of  the  transition  from  one  stage  to  the  next  may 
affect  the  results  of  later  stages  of  calculation. 
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CALCULATIONAL  MATRIX  TAMPED  150  KT  EXPLOSIONS 
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Based  on  extensive  experience,  a multistage  operational 
sequence  was  adopted  for  executing  the  original  matrix  of  ground- 
motion  problems.  It  consisted  of  the  following  seven  primary 
stages  of  calculation: 

1.  Calculation  of  the  motion  from  bang  time  until  the 
leading  shock  produced  by  the  explosion  propagates 

a distance  almost  as  great  as  the  depth  of  burial  or 
to  a radius  several  times  the  final  cavity  radius. 
During  this  stage  the  field  is  assumed  to  be  spher- 
ically symmetric  and  the  calculation  is  performed 
using  the  AFTON  1 code. 

2.  Generation  of  initial  conditions  for  the  AFTON  2A 
code  from  the  AFTON  1 field  at  the  close  of  Stage  1, 
followed  directly  by  the  calculation  (using  AFTON  2A) 
of  axisymmetric  motion  during  the  early  period  of 
crater  formation. 

3.  Redistribution  of  mesh  points  to  form  a rectangular 
array  when  cracked  rock  about  the  burst  point  has 
become  so  porous  that  mixing  of  cavity  gases  with 
rock  is  more  probable  than  separation  of  the  two  by 
a sharp  boundary. 

4.  Resumption  of  the  calculation  in  the  new  finite- 
difference  mesh  with  the  addition  of  cells  to 
accommodate  the  expanding  field  of  motion.  Calcula- 
tion continues  until  fine  definition  of  the  field 

in  the  vicinity  of  the  crater  is  no  longer  needed. 

5.  Removal  of  alternate  vertical  and  horizontal  rows 
of  mesh  points  to  form  a coarser  finite-difference 
mesh.  The  generalized-coordinate  capability  of  the 
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AFTON  2A  code  is  used  in  Stage  5 so  that  mass, 
momentum,  kinetic  energy,  and  internal  energy  are 
all  exactly  conserved. 

6.  Further  addition  of  cells  to  the  finite-difference 
mesh  so  that  the  mesh  continues  to  enclose  the  grow- 
ing region  of  disturbed  material. 

7.  Deletion  of  cells  from  the  outer  regions  of  the 
mesh  to  contract  the  domain  of  computation.  The 
rate  of  contraction  is  chosen  to  ensure  that  the 
field  is  calculated  to  a minimum  time  of  2.5  sec  of 
motion  on  a certain  hemisphere  without  contamina- 
tion by  signals  from  the  non-physical  contracting 
mesh  boundary.  This  hemisphere,  whose  radius 
defines  a so-called  "magic  circle"  in  the  azimuthal 
plane  of  calculation,  is  used  as  a source  surface 
for  subsequent  calculations  of  teleseismic  radia- 
tion outside  the  region  of  nonlinear  material 
behavior . 

C3.  Difficulties  Encountered  while  Processing  the  Results  of 

the  Calculations 

During  calculation,  selected  data  were  stored  as  functions 
of  time  to  be  processed  after  completing  the  source  calcula- 
tion. The  data  included  a description  of  material  ejected 
through  the  ground  surface  as  well  as  field  variables  at  one- 
degree  intervals  along  two  magic  circles  — quarter-circles, 
actually,  centered  at  ground  zero  in  the  plane  of  motion.  The 
data  on  ejecta  were  used  to  predict  crater  size  (Appendix  D) , 
unless  burial  depth  was  great  enough  to  prevent  ejecta  produc- 
tion and  crater  formation.  The  magic-circle  data  were  at  first 

3 

obtained  for  use  by  Systems,  Science  and  Software  (S  ) in 


computing  body  and  surface  wave  magnitudes,  m,  and  M , for 

these  calculated  fields  of  motion.  Late  in  the  program,  Pacific 

Sierra  Research  Corporation  (PSRC)  made  use  of  the  same  data  to 

provide  an  alternative  set  of  values  of  Mg . While  pursuing 

calculation  of  surface-wave  magnitudes,  PSRC  found  that  a net 

vertical  impulse  had  been  delivered  to  the  hemispherical  region 

bounded  by  the  ground  surface  and  a quarter-circle  in  the  plane 

of  motion.  This  impulse  made  a significant  contribution  to  M . 

s 

Since  no  sources  were  located  outside  the  region  and  the  ver- 
tical momentum  for  the  material  within  the  region  was  almost 
zero,  the  presence  of  the  non-zero  vertical  impulse  found  by 
PSRC  suggested  that  our  computational  medium,  prior  to  explosive 
disturbance,  was  not  in  exact  static  equilibrium.  The  net 
impulse  delivered  to  the  hemispherical  region,  which  contains  all 
sources,  can  be  calculated  by  integrating  the  time-dependent 
surface  tractions  and  body  forces.  In  the  presence  of  a gravita- 
tional field  the  body  force  is  the  weight  of  the  material  occupy- 
ing the  region,  i.e.,  the  vertical  force  W^m^g.  The  vertical 

force  due  to  the  atmospheric  pressure  applied  to  the  top  surface 

2 

of  the  closed  hemisphere  is  given  by  At=P  rrR^  where  the  sub- 
script t indicates  the  time  at  which  the  force  is  applied  to  the 

surface.  Finally,  the  integral  of  the  surface  traction  over  the 

t f t 

curved  surface  is  given  by  S^=  nj  ds  > where  n^  (j  - 1,  2 or  3) 

is  the  outward  normal.  The  atmospheric  force  and  the  gravita- 
tional force  have  only  vertical  components.  As  a consequence 
of  symmetry,  the  surface-traction  integral  contributes  only  a 
vertical  component  S^.  Therefore  vertical  impulse  is  given  as  a 
function  of  time  by 


(Cl) 
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The  integral  vanishes  at  t=0.  Since  the  material  inside  the 
region  is  initially  at  rest,  Iq=0.  Therefore  conservation  of 
momentum  implies  the  integral  is  zero  at  t=®  when  the  material 
is  once  again  at  rest.  (This  is  very  nearly  the  case  at  the 
completion  of  the  calculation  when  t = 2.5  seconds.) 

To  allow  credible  calculation  of  surface  wave  magnitudes 

and  other  measures  of  the  effect  of  these  bursts,  the  source  of 

the  discrepancy  between  the  impulse  calculated  by  PSRC  and  the 

momentum  calculated  by  ATI  had  to  be  found  and  corrected.  A 

test  problem  was  run  to  determine  if  the  equations  in  the  AFTON 

code  did  indeed  conserve  mass,  momentum  and  energy.  The  output 

of  that  calculation  showed  that  the  conservation  laws  were  being 

obeyed.  Among  the  quantities  calculated  for  these  conservation 

checks  were  the  forces  acting  along  the  thermodynamic  cell 

boundaries  (which  make  up  the  surface  of  the  region)  and  the 

velocities  of  the  mesh  points.  These  stress  and  velocity  data 

are  consistent  with  the  discrete  conservation  relations  in  the 

AFTON  code  and  are  the  quantities  to  be  used  in  calculating  M . 

s 

The  data  originally  provided  to  PSRC  represented  stress  at 
specific  target  points,  not  the  average  stress  acting  on  a cell 
boundary.  It  was  therefore  inappropriate  for  calculating 
impulse.  Compatible  thermodynamic  cell  stress  and  mesh  point 
velocity  data  were  therefore  obtained  during  the  running  of 
subsequent  problems. 

Since  the  importance  of  vertical  impulse  in  the  calcula- 
tion of  M had  been  demonstrated,  details  of  the  ground  motion 
s 

calculation  procedure  were  reviewed  to  identify  and  remove  as 
many  sources  of  numerical  inaccuracy  affecting  vertical  impulse 
as  possible.  Among  the  numerical  changes  required  were  1)  a 
redefinition  of  the  midpoint  of  a thermodynamic  cell  and  2) 
revision  of  the  procedure  for  balancing  the  gravitational  field. 
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The  midpoint  definition  was  changed  from  one  which  attempted  to 
balance  the  mass  in  the  four  quarter  cells  to  one  which  gave  a 
simple  relationship  between  the  volume  and  area  of  rectangular 
cells.  As  long  as  the  undisturbed  zones  are  rectangular, 
calculation  of  the  geostatic  stress  field  can  be  based  on  a 

l 

simple  algorithm.  A small  correction  is  made  to  the  original 
definition  of  the  gravitational  stress  field. 

The  largest  numerical  error  in  balancing  the  gravita- 
tional forces  occurs  when  the  cells  are  arbitrary  quadrilaterals. 
Although  these  errors  are  less  than  0.2%  of  the  gravitational 
acceleration,  the  errors  introduce  displacements  of  the  same 
order  as  the  correct  permanent  displacements  (Section  1.3). 

The  acceleration  errors  in  non-rectangular  cells  rise  from  the 
unrelated  procedures  for  calculating  forces  and  masses  from  cell 
shape.  As  cell  shape  departs  further  and  further  from  the 
desired  rectangular  form,  midpoint  depths  in  adjacent  cells  may 
begin  to  differ,  introducing  additional  sources  of  error.  Signif- 
icant numbers  of  non-rectangular  cells  are  used  in  the  initial 
AFTON  2A  mesh  (a  discrete  vapor-solid  interface  must  be  carried 
during  this  stage  of  calculation) . A method  of  balancing  the 
vertical  forces  was  therefore  necessary. 

Cl 

As  described  in  the  AFTON  users  manual,  the  equation 
for  the  change  in  momentum  M during  time  step  Atn  = tn  - tn  ^ is 

M"  - Mn_1  = At  |(Fla  + Flb  + Flc  + Fld)  + body  forces|  (C2) 

where  a subscripted  number  indicates  a mesh  point  and  the  sub- 
scripted letters  indicate  the  thermodynamic  cells  having  that 
meshpoint  in  common.  Starting  with  zone  'a'  above  and  to  the 
left  of  the  meshpoint,  the  remaining  zones  are  encountered  in 
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clockwise  order.  By  starting  with  the  cell  adjacent  to  the  axis 
of  symmetry  and  just  below  the  ground  surface  and  proceeding 
horizontally  to  the  mesh  boundary  and  then  returning  to  the  axis 
of  symmetry  for  the  next  row  of  cells,  etc.,  the  force  and  mass 
in  zone  'c'  can  be  calculated  from  the  values  already  calculated 
for  three  other  cells  according  to  the  relation: 
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The  stress,  strain  (P,u)  state  in  the  statically  balanced 
gravational  field  is  then  obtained  from  an  iterative  solution  of 

F*  = g PQ(1  + n)Vl(;  + P(u)A^c  (C4) 

where  p is  the  reference  density  of  the  material,  V,  is  the 
o J * lc 

volume  of  the  quarter  cell  and  A^c  is  the  vertical  component  of 
the  area  of  the  quarter  cell.  The  associated  horizontal  forces 
are  generally  not  balanced  and  a net  horizontal  acceleration 
results.  To  correct  this,  a pseudo-acceleration  is  introduced. 
The  pseudo-acceleration,  varying  from  cell  to  cell,  is  deter- 
mined by  balancing  the  horizontal  momentum  equation;  it  appears 
as  a constant  in  the  subsequent  calculation  of  axisymmetric , 
source-induced  ground  motion. 

To  determine  if  the  initial  field  obtained  in  this  way 
is  actually  in  static  equilibrium,  a calculation  was  made  in 
which  there  were  no  source  terms.  The  average  velocities, 

£m,v,  . 

V.  * -=^ — ■*-,  in  the  field  after  100  timesteps  of  computation 
i 2^m . 
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were  less  than  10  times  the  gravitational  velocity,  V =gt  in 
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the  vertical  direction  and  less  than  10  cm/sec  in  the  horizon- 
tal direction.  The  velocities  were  randomly  oriented  and  no 
systematic  effect  of  the  gravitational  field  could  be  found. 

C4 . Numerical  Errors  Caused  by  Calculational  Technique 

In  the  original  calculations  the  first  change  in  zone 
configuration  was  a rezone  which  gave  all  cells  the  same 
horizontal  dimensions  and,  except  for  the  zones  just  below  the 
ground  surface,  the  same  vertical  dimensions.  To  achieve  this 
configuration  change,  most  lines  of  mesh-points  were  moved  large 
distances,  by  a series  of  finite  displacements  at  zero  timestep. 
This  grid  motion  (and  the  accompanying  transport)  tended  to  smear 
the  signal:  after  the  rezone  pulse  shape  had  been  altered. 

Pulse  amplitudes  were  reduced,  high  frequency  content  was  atten- 
uated, and  arrival  times  of  the  direct  and  reflected  waves  were 
altered.  General  pulse  broadening  occurred.  On  a cycle-by-cycle 
basis,  the  mass,  momentum,  kinetic  energy  and  internal  energy 
were  conserved  locally  during  this  rezone,  whence  they  were  also 
conserved  overall.  Material  transport  during  grid  line  movement 
did,  however,  cause  a net  redistribution  of  energy  and  momentum 
in  the  field. 

Examination  of  the  original  computational  mesh,  the  mesh 
desired  after  rezone  and  the  grid  line  movement  necessary  to 
accomplish  the  mesh  configuration  change  led  to  two  modifica- 
tions. First,  the  graduated  spacing  of  grid  lines  and  the  non- 
rectangular  cell  shape  in  cells  away  from  the  burst  point  in 
the  original  mesh  were  replaced  with  uniform  rectangular  spacing. 
Second,  the  rezone  procedure  was  modified.  (This  was  strongly 
influenced  by  the  modified  original  mesh.)  The  new  technique 
calls  for  the  removal  of  about  1/3  of  the  grid  lines  in  each 
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direction.  Since  the  modified  original  grid  is  made  up  of  cells 
which  are  approximately  square,  this  technique  implies  that 
post-dezone  cell  boundaries  will  have  moved  at  most  a cell  width 
and  usually  only  half  a cell  width.  In  the  case  where  graduated 
spacing  was  used,  grid  lines  moved  several  cell  widths.  As 
before,  mass,  momentum,  internal  energy  and  kinetic  energy  are 
conserved  both  for  the  whole  grid  and  locally  on  a cycle-by- 
cycle basis.  However,  since  material  transport  for  any  cell 
involves  only  a small  region  around  its  original  location  in  the 
field,  redistribution  of  momentum  and  energy  is  limited  to  local 
regions  — which  improves  the  accuracy  of  numerical  solutions. 

C5.  Results  of  the  Changes 

Test  problems  were  run  to  check  the  effect  of  each 
change.  After  all  changes  were  completed,  5 of  the  original 
problems  were  recalculated  and  4 additional  calculations  were 
made.  Comparison  of  body  wave  magnitudes  obtained  from  the 
original  and  revised  computation  procedures  (Table  CT2)  showed 
only  small  changes  except  for  Problem  12,  Wet  Sandstone  @ DOB=531.3m 
The  large  change  observed  in  Problem  12  is  due  to  the  discontiu- 
ity  in  the  seismometer  displacement  definition  of  m^  (Sec- 
tion 3.3).  When  the  epoch  of  seismometer  responses  determining 
m^  contains  features  such  as  an  emerging  local  extremum  which 
may  cause  a jump  in  the  period,  a sudden  change  in  peak-to-peak 
amplitude  or  the  switching  of  the  maximum  from  one  excursion  to 
another,  small  changes  in  the  driving  function  can  cause  jumps 
in  m^.  For  the  problems  treated  here,  seismometer  response  is 
contaminated  by  the  inclusion  of  reflected  S-wave  signals.  When 
there  is  little  separation  between  the  P-  and  S-wave  contribu- 
tions to  the  ground  motion,  seismometer  response  and  body-wave 
magnitude  can  be  very  sensitive  to  small  changes  in  the 
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Table  C2 


t No  cracking  or  gravitational  effects  in  calculation 
**No  gravitational  effects  in  calculation 
-M-No  cracking  in  calculation 


computational  procedure.  Because  of  the  differences  in 
wavespeeds,  separation  of  P-  and  S-wave  signals  can  be  obtained 
by  observing  ground  response  at  greater  ranges.  By  running  the 
problems  to  later  times  and  monitoring  particle  velocity  at 
those  distant  stations,  the  effect  of  S-wave  contamination  on 

m,  can  be  eliminated. 

b 

Body-wave  magnitudes  were  computed  from  calculated  near- 
field source  motions  observed  at  slant  ranges  between  0.7  and 
5.8  km  for  six  bursts  in  wet  sandstone,  Problems  10,  12,  26, 

110,  jl12  and  126.  Values  of  m^  from  observation  stations  along 
the  lines  for  0=75°  and  9=85°  (see  Figures  C-l(a))  were  averaged 
and  plotted  as  functions  of  slant  range  R^  to  illustrate  the 
effect  of  wave  separation.  For  those  cases  in  which  the  effects 
of  gravity  and  cracking  w^-e  included  (Problems  10,  12  and  26), 
separation  of  P-  and  S-wave  response  in  the  near-fieid  is  effec- 
tively completed  at  observational  slant  ranges  of  about  3 km, 
Figure  C-l(b).  The  discontinuity  in  associated  with  sudden 

changes  in  pulse  shape,  period,  or  time  of  occurrence  of  max- 
imum peak-to-peak  seismometer  displacement  is  illustrated  by  the 
dashed  curves  AB'C  and  ABC  for  Problem  12.  At  some  range  between 
0.7  and  1.78  km  for  0=75°  (between  1.78  and  2.8  km  for  0=85°), 
the  definition  of  mb  is  transferred  discontinuous ly  (Section  3.3) 
from  one  identifiable  feature  to  another. 

For  the  cases  where  the  effects  of  gravity  and  cracking 
were  omitted  (Problems  110,  112,  and  126;  Figure  C-l(c)),  the 
close-in  variation  of  magnitude  with  slant  range  is  DOB  depen- 
dent. For  shallow  burial,  Problem  110  @ 207.2m,  mb  shows  an 
initial  decrease  with  slant  range.  For  greater  depths  of  burial, 
an  initial  increase  is  observed.  In  all  three  cases  little  varia- 
tion in  m^,  and  especially  in  relative  values  of  m^,  occurs  for 
R >3km.  Magnitudes  obtained  from  source  pulses  observed  at 
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nominal  slant  ranges  of  5 km  are  compared  with  those  for  source 
pulses  from  nominal  slant  ranges  of  1 km  in  Table  C2.  On  the 
basis  of  these  six  problems,  it  appears  that  for  wet  sandstone 
the  effect  of  shear  wave  contamination  on  m^  can  be  avoided  by 
using  as  teleseismic  sources,  velocity  pulses  from  slant  ranges 
greater  than  3 km. 

C6.  Revised  Calculational  Technique 

After  finding  a way  to  achieve  exact  static  equilibrium 
in  discrete  fields,  and  identifying  the  effects  of  several  varia- 
tions in  procedure,  a revised  calculation  technique  was  adopted. 
The  revised  technique  retains  the  seven  basic  stages  of  the 
original  procedure  (Section  C2)  but  utilizes  modified  routines 
for  balancing  the  gravitational  field  and  modified  computational 
meshes  and  dezone  procedures.  Consistent  stress  and  velocity 
data  are  obtained  for  use  in  computing  surface-wave  amplitude, 
and  stations  used  to  obtain  teleseismic  source  data  (velocity 
pulses)  for  body-wave  magnitude  calculation  are  moved  to  ranges 
where  P-  and  S-wave  contributions  are  separated  in  time. 

Stage  1,  calculation  of  the  axisymmetric  phase  of  early 
motion,  is  unaffected  by  the  revisions.  Redesign  of  the  grid 
used  to  initiate  the  AFTON  2A  calculation  of  axisymmetric  motion 
affects  Stage  2.  In  the  revised  procedure,  outside  a small 
region  which  includes  the  burst  point  and  the  zone  of  strong 
deformation  associated  with  cavity  expansion,  the  computational 
grid  (Grid  I)  is  composed  of  nearly  uniform  rectangular  cells. 
Originally,  graduated  spacing  had  been  used  with  smaller  zones 
near  the  burst  point  to  resolve  close-in  phenomenology  and 
larger  cells  at  larger  ranges  where  slower  variations  in  space 
and  time  are  observed.  Before  the  initial  conditions  are 
inserted  and  the  axisymmetric  calculation  begun,  the  precise 
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stresses  required  for  static  equilibrium  of  the  host  medium  in 
a gravitational  field  are  calculated  for  Grid  I in  its  entirety 
(Section  C3).  The  initial  grid  is  then  activated,  initial 
conditions  introduced,  and  calculation  begun.  As  the  disturbed 
region  expands,  additional  cells,  already  carrying  the  stresses 
needed  to  balance  the  gravitational  body  forces,  are  activated. 
Calculation  proceeds,  with  the  mesh  points  kept  Lagrangian 
except  near  ground  zero  where  ejecta  production  is  being 
monitored,  until  the  activation  of  a cell  at  a depth  of  1200m. 
(For  the  largest  depth  of  burial  - 743.9m  — this  limiting  depth 
was  increased  to  1750m).  Stage  3,  a dezone,  is  then  initiated. 

Dezoning  is  accomplished  in  two  steps:  1)  Grid  lines 
to  be  removed  are  moved  in  non-Lagrangian  fashion,  at  a nearly 
vanishing  timestep,  to  positions  almost  coincident  with  grid 
lines  to  be  retained.  2)  Grid  lines  are  removed  and  material 
in  the  cells  about  to  vanish  is  'mixed*  with  material  in  the 
cells  to  be  retained.  First,  horizontal  lines  are  removed, 
then,  by  repeating  the  two  steps,  vertical  lines  are  removed. 

In  these  calculations  two  kinds  of  dezone  are  employed:  selec- 
tive dezones  and  alternate  dezones.  Stage  3 is  normally  a 
selective  dezone.  In  this  case  all  grid  lines  move  toward  grid 
line  positions  desired  after  the  dezone:  lines  to  be  retained 
move  to  exactly  those  positions;  lines  to  be  removed  — about  1/3 
of  the  lines  in  either  direction  — approach  those  positions. 
Except  for  the  ground  surface  all  new  grid  lines  are  straight. 

In  the  region  surrounding  the  burst  point  and  extending  upward 
through  the  crater  region,  movement  and  straightening  of  grid 
lines  during  the  selective  dezone  causes  material  transport 
across  cell  boundaries  and  mixing  of  material  in  very  different 
states.  Elsewhere,  the  effects  of  grid  line  movement  are  less 
severe.  Revision  of  the  grid  to  the  form  of  Grid  I and  the  use 


of  selective  dezoning  at  Stage  3 reduces  the  disruptive  influence 
by  minimizing  the  translation  of  grid  lines  and  localizing  the 
transport  and  mixing.  For  these  problems  grid-line  motion  is, 
for  the  most  part,  limited  to  one  cell  width.  After  a selec- 
tive dezone,  zone  volume  has  increased  to  about  2\  times  the 
original  cell  volume.  Mass,  momentum,  and  energy  are  conserved 
and  the  time  step  is  automatically  increased. 

In  alternate  dezones  (Stage  5) , every  other  vertical 
grid  line  and  every  other  horizontal  grid  line  is  removed  to 
produce  new  cells  with  volumes  about  4 times  the  original  cell 
volume.  The  alternate  dezone  requires  the  same  two-phase,  two- 
step  procedure  of  moving,  then  removing  lines.  However,  the 
lines  to  be  retained  are  not  moved. 


Following  a dezone,  a new  grid  is  constructed  and 
gravitational  body  forces  balanced.  During  Stage  3,  the 
active  portion  of  dezone  Grid  I is  used  as  the  core  for 
constructing  Grid  II.  Similarly  during  Stage  5 the  active 
part  of  dezoned  Grid  II  is  used  as  the  core  for  Grid  III.  For 
Grid  II,  nearly  uniform  rectangular  cells  are  added  along  the 
bottom  and  outer  radius  of  the  core  grid.  The  gravitational 
body  force  is  balanced,  an  initial  grid  carrying  initial  condi 
tions  from  dezoned  Grid  I is  activated,  and  the  calculation 
resumed.  Additional  cells  are  activated  (Stage  4)  and  calcula 
tion  continues  until  the  next  dezone  is  required.  The  normal 
sequence  calls  for  a selective  dezone  and  a later  alternate 
dezone,  but  under  special  circumstances  the  type  and  number  of 
dezones  may  change.  For  example,  for  the  shallowest  burst  in 
granite  (DOB-159m),  three  dezones  were  used  in  the  sequence: 
alternate,  alternate,  selective,  to  allow  proper  resolution  of 
events  in  the  crater  region. 


The  monitoring  of  ejecta  data  at  the  ground  surface 
begins  during  Stage  2 and  continues  through  problem  termination. 
However  (see  Appendix  D)  to  avoid  perturbations  in  ejecta  mass 
caused  by  the  dezones,  ejecta  data  prior  to  the  second  dezone 
and  data  from  the  field  of  motion  within  the  developing  crater 
zone  at  the  close  of  Stage  2,  form  the  primary  basis  for  describ 
ing  the  completed  crater.  Stations  for  monitoring  velocity  and 
compatible  stress  data  on  the  hemispherical  surface  outside  the 
region  of  inelastic  response  are  normally  established  for  sand- 
stone after  the  first  dezone.  Since  further  dezoning  will  occur 
before  problem  termination,  care  must  be  taken  in  defining  the 
hemisphere  (i.e.,  the  quarter  circle  in  the  mesh  or  r,z-plane) 
so  that  it  persists  as  a continuously  identifiable  feature 
through  the  dezone.  For  bursts  in  granite  the  hemispherical 
surface  for  monitoring  stress  and  velocity  data  is  not  estab- 
lished until  Stage  5 (following  the  second  or,  for  the  shallow- 
est burst,  third  dezone). 

Stage  5,  the  last  dezone,  is  initiated  for  bursts  in 
sandstone  by  activation  of  cells  at  the  boundary  of  Grid  II  and 
in  granite  by  signal  arrival  on  a hemisphere  of  radius  2100m. 

The  monitoring  stations  on  the  hemisphere  in  sandstone  are 
redefined  in  terms  of  Grid  III  which  is  built  around  a core 
from  dezoned  Grid  II.  In  granite  the  monitoring  stations  on 
the  hemisphere  are  finally  defined.  Calculation  then  proceeds 
with  the  activation  of  additional  cells.  Stage  6,  until  tho 
limit  of  Grid  III  is  reached.  At  that  point,  Stage  7 begins 
with  deletion  of  grid  lines  at  the  outer  boundary,  propagating 
inward  at  a rate  faster  than  the  local  longitudinal  wavespeed, 
to  prevent  spurious  signals  reflected  from  the  artificial  mesh 
boundary  from  entering  the  calculation.  The  dropping  of 
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vertical  and  horizontal  lines  is  handled  independently. 

Stage  7 continues  until  problem  termination  time,  t-2.5  seconds, 
when  the  active  grid  has  shrunk  almost  down  to  the  hemispherical 
monitoring  surface. 
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Dl.  Phenomenology 

Near  surface  nuclear  explosions  produce  craters  through 
the  action  of  three  primary  mechanisms:  ejecta  production, 
plastic  flow,  and  compaction.  Crater  size  and  the  relative  im- 
portance of  each  of  these  mechanisms  changes  with  depth  of 
burst  (DOB)  and  to  a lesser  degree  with  site  medium.  Below  some 
DOB,  ejecta  production  ceases.  For  this  investigation  bursts  at 
or  below  that  depth  are  not  considered  to  produce  craters. 

(Three  of  the  original  twelve  problems  had  burst  points  below 
that  depth;  no  crater  data  are  given  for  those  cases.) 

A typical  crater  is  shown  in  cross  section  in  Figure  Dl. 
The  major  features  shown  are:  the  apparent  crater  (the  visible 
limit  of  the  "hole  in  the  ground") , the  true  crater  (the  limit 
of  dissociation  of  the  natural  material),  fallback  (broken  ma- 
terial filling  the  volume  between  the  apparent  and  true  craters), 
the  rupture  zone  in  which  cracking  (but  not  dissociation)  occurs, 
a plastic  zone  where  less  severe  inelastic  deformations  occur, 
and  finally,  an  elastic  zone.  For  the  DOB's  considered  here, 
ejecta  production  is  the  dominant  crater  forming  mechanism. 

Ejecta  comes  almost  exclusively  from  material  within  a cone 
whose  apex  lies  at  the  burst  point  and  whose  surface  extends 
upward  to  the  ground  surface.  In  competent  rock  or  saturated 
media  the  cone  half  angle  is  approximately  45°,  i.e.,  the  limits 
of  the  ejecta  cone  are  approximately  those  of  a cone  of  shear 
failure  reaching  the  surface.  Sketches  of  true  craters  from 
buried  explosions0*  support  this  correlation. 

Outside  the  ejecta  cone,  material  will  flow  away  from  the 
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burst  point  during  the  early  phases  of  crater  formation  when 
stress  levels  are  high  enough  to  cause  shear  failure.  Later, 
when  stress  levels  have  fallen,  subsequent  deformation  will 
be  elastic.  At  this  point  the  resistance  of  the  material  to 
elastic  shear-strain  makes  further  displacement  insignificant 
on  the  scale  of  crater  size.  Thus,  the  passage  of  material 
near  the  crater  region  from  states  of  shear  failure  to  elastic 
states  resembles  the  freezing  of  a fluid  which  had  been  under- 
going large  distortions.  "Freezing"  of  the  material  as  plastic 
flow  ceases,  effectively  halts  further  growth  of  the  true  crater 
at  a time  that  varies  from  a few  tenths  of  a second  to  a few 
seconds  for  the  bursts  of  interest  here,  depending  on  the  prop- 
erties of  the  medium  and  on  the  energy  released  in  the  burst. 


Compaction,  i.e.,  permanent  volume  change  under  purely 
hydrostatic  loading,  is  of  no  consequence  for  buried  explosions 
in  NTS  granite,  Soviet  granite,  and  wet  sandstone.  In  these 
materials  saturation  or  the  natural  lack  of  porosity  led  to 
material  models  with  no  air  filled  void  space,  and  consequently 
no  compaction.  For  the  dry  sandstones  (dry  sandstone  and  weak 
dry  sandstone)  significant  compaction  did  occur  in  the  close-in 
region.  Its  effects  due  either  to  the  rapid  attenuation  of  the 
outgoing  signal  (energy  dissipation)  or  to  the  permanent  volume 
reduction  (compaction) , is  automatically  accounted  for  through 
the  density  of  each  cell. 


Two  other  mechanisms  play  a role  in  determining  crater 
size:  fallback  and  bulking  of  material  within  the  rupture  zone 
In  the  latter  phases  of  crater  formation  and  for  as  much  as  a 
few  tens  of  seconds  thereafter,  ejecta  will  rain  down  on  the 
ground  surface  both  inside  and  outside  the  crater.  Within  the 
true  crater  this  material,  fallback,  produces  so  thick  a layer 


of  rubble  that  the  depth  of  the  apparent  crater  left  by  the 
burst  is  generally  no  more  than  half  that  of  the  true  crater. 

In  the  rupture  zone  deformation  of  the  medium  is  severe  enough 
to  cause  cracking.  This  cracking,  a local  phenomenon  which 
expands  the  material  without  causing  dissociation,  modifies 
crater  volume  by  swelling  the  the  immediate  surroundings  of 
the  crater. 

D2.  Determination  of  Crater  Size 

Four  phenomena  were  considered  in  determining  crater 
size:  1)  ejecta  production,  2)  plastic  volume  change,  3)  fall- 
back, and  4)  bulking.  Contributions  to  crater  volume  were  com- 
puted for  each  phenomenon  and  combined  to  obtain  estimates  of 
crater  volume.  Crater  radius  was  used  to  characterize  crater 
size.  By  assuming  a standard  shape  a simple  relationship  be- 
tween crater  radius  and  crater  volume  was  determined.  On  the 
basis  of  crater  profiles  for  several  Soviet  peaceful  nuclear 
explosive  (PNE)  events0^",  an  inverted  truncated  conical  crater 
shape  was  adopted.  The  data  indicate  that  the  crater  bottom  is 
nearly  flat  out  to  a fraction,  «(“0.4),  of  the  apparent  crater 
radius  Rc>  The  walls  of  the  crater  slope  upward  at  an  angle 
0(«3O°).  The  volume  of  the  assumed  standard  crater  is 

3 

V 3 

V£  - — ^ — (1-or  )tan0 

For  a«0.4  and  0-30°,  the  volume  is 

v „ Q.936tt  r3 
c 3 v/3  c 

The  crater  radius  and  depth  Dc  are  given  by 
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3 v r * 

R - == = 1.209  V J ; D = (l-er)R  Tan8  = 0.3464  R 

C |Tr(l-0f3)taneJ  c c c c 

Each  contribution  to  crater  volume  is  obtained  as  a func- 
tion of  range.  By  plotting  the  estimated  crater  volumes  (obtained 
by  summing  contributions  from  a prescribed  set  of  mechanisms)  as 
functions  of  range  on  the  crater  volume,  crater  radius  plot  as 
shown  in  Figure  D2,  estimates  of  crater  radius  are  obtained  from 
the  intersections  (Points  1,2, 3, 4 and  5). 

D3.  Details  of  Crater  Size  Calculation 

To  account  for  effects  of  the  four  mechanisms  contribut- 
ing to  cratering,  problem  output  was  monitored  to  determine  the 
mass,  speed,  direction  and  time  of  departure  of  ejecta  particles 
for  subsequent  calculation  of  their  ballistic  trajectories;  the 
cells  in  which  shear  failure  occurs;  and  the  cells  in  which 
cracking  occurs.  Using  this  data  and  data  normally  used  in  ad- 
vancing the  ground-motion  calculation,  the  crater  volume  contri- 
butions from  each  mechanism  were  calculated  as  described  below. 

D3.1  Calculations  of  Ejecta 

To  calculate  the  amount  of  material  ejected  from  the 
ground  in  the  original  twelve  calculations  of  axisymmetric  mo- 
tion, a special  boundary  condition  was  applied  at  the  ground 
surface.  The  upward  jump-off  motion  at  the  ground  surface  could 
be  of  two  types  1)  ejecta  production  and  2)  surface  spall.  For 
ejecta  production,  the  vertical  velocity  assumes  a very  large 
value  and  the  material  is  broken  up  into  many  pieces  of  varying 
size.  The  exit  velocity  is  sufficient  to  carry  this  material 
large  distances  from  its  point  of  origin  before  it  returns  to 

the  surface  after  a ballistic  flight.  For  surface  spall,  the 

0 
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Figure  D2:  Comparison  of  Crater  Volumes  Estimated  from  Computed  Ground 
Motion  with  Volumes  of  Truncated  Cone  Craters.  Intersection 
of  Curve  with  "truncated  cone"  curve  determines  crater  radius 


vertical  velocities  are  relatively  small.  The  material  cracks 
but  remains  fairly  competent;  it  returns  to  the  original  ground 
surface  level  very  near  its  place  of  origin.  The  total  movement 
of  spall  material  is  only  a few  meters. 

To  differentiate  between  ejecta  and  spall  material,  a 
set  of  screening  criteria  were  developed.  To  be  counted  as 
ejecta  material  just  below  the  ground  surface  had  to  fulfill 
three  requirements:  1)  the  vertical  component  of  particle  veloc- 
ity is  positive,  2)  the  current  position  of  the  ground  surface 
is  above  the  original  ground-surface  level,  and  3)  the  material 
just  below  the  ground  surface  has  expanded  so  that  its  density 
is  less  than  90  percent  of  the  original  density,  p.  When  all 
three  of  these  conditions  are  satisfied,  the  elevation  of  the 
ground  surface  becomes  fixed  in  space  and  material  is  transported 
through  it.  If  any  of  the  three  conditions  are  not  met,  the 
ground  surface  moves  in  a Lagrangian  manner  and  no  material  is 
ejected.  In  nine  of  the  twelve  calculations,  material  was  ejected 
using  the  conditions  described  above.  The  three  calculations  that 
produced  no  ejecta  were  the  bursts  in  dry  sandstone  @ 253.0m 
(Problem  6)  and  in  wet  sandstone  @ 253.0  and  531.3m  (Problems  11 
and  12  respectively) . 

Although  monitoring  ejecta  theoretically  provides  an  ac- 
curate account  of  total  ejecta  mass,  the  dezones  necessary  to 
run  the  calculation  introduce  strong  perturbations,  leaving  the 
accuracy  of  the  total  mass  in  doubt.  To  avoid  this  problem, 
ejecta  mass  was  monitored  in  the  manner  just  described  through 
the  time  when  the  first  dezone  was  made.  At  that  time,  the  field 
variables  were  used  to  determine  which  material  in  the  ground 
would  be  ejected  — the  criterion  now  being  that  particles  below 
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the  surface  would  have  to  have  sufficient  vertical  velocity  to 
clear  the  ground  surface  under  the  effect  of  gravity.  Although 
the  density  of  this  material  in  the  ground  was  higher  than  that 

: 

needed  to  form  ejecta  at  the  ground  surface,  it  was  assumed 
that  by  the  time  the  particle  reached  the  ground  surface  its 
density  would  satisfy  the  ejecta  density  criterion.  This  mod- 
ification tends  to  overestimate  ejecta  mass. 

To  determine  the  mass  (volume)  of  material  ejected  as  a 

function  of  range,  the  ground  surface  is  divided  into  concentric 

, A 

rings  of  equal  width,  (5(W(KT)/150)^) . The  total  mass  exiting 
(or  which  will  exit  after  the  first  dezone)  through  each  ring, 
mej , is  divided  by  the  initial  density,  pq,  to  get  the  volume 
ejected  through  that  ring.  Summing  the  volumes  ejected  through 
each  ring  from  the  center  out  to  the  range,  r.,  then  gives 
ejecta  volume  as  a function  of  range. 

D3.2  Calculation  of  Plastic  Volume  Change 

The  next  phenomenon  considered  is  volume  change  during 
plastic  flow.  To  determine  the  volume  change  due  to  plastic 

I flow,  V the  mass  of  material  in  an  annulus  of  material  ex- 

< * p£’ 

tending  from  the  surface  down  to  a depth,  Yq,  where  no  signals 
had  arrived,  was  determined  from  the  density  profile  at  a ref- 
erence time.  The  height  that  material  would  have  if  it  were  at 
the  initial  density  was  determined  by  dividing  the  mass  by  the 
initial  density,  p , and  the  area.  This  height  was  subtracted 
from  Yq  to  give  the  height,  hpg . , in  that  annulus  occupied  by 
the  material  ejected  and  the  material  removed  due  to  plastic 
flow.  To  correct  for.  ejecta  that  leaves  the  annulus  before  the 
reference  time,  the  height  of  an  annulus  of  the  already  ejected 

% material  was  found  and  subtracted  from  h . . The  corrected  con- 
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tribution  of  plastic  flow  to  crater  volume  for  an  annulus  is  then 


X 


The  values  for  each  annulus  are  then  summed  from  ground  zero  out 
to  radius  r^  to  obtain  plastic  volume  change  as  a function  of 
range. 

D3.3  Fallback 


In  conjunction  with  the  monitoring  done  to  determine  ejec- 
ta mass  as  a function  of  time  and  position,  sufficient  data  were 
also  obtained  (speed,  direction,  position  and  time  of  departure) 
to  allow  calculation  of  the  ballistic  flight  of  each  ejecta  mass. 
Using  the  ejecta  starting  positions  and  ballistic  trajectory  data, 
the  mass  of  fallback  material  accumulating  in  each  ring  is  deter- 
mined. Since  the  fallback  accumulates  as  a loose  blanket  of  rub- 
ble, its  density  pf«fp  is  less  than  the  density  of  the  original 
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undisturbed  material  (f<l) . For  typical  hard  rock  , the  density 

of  fallback  is  0.78pQ.  The  volume  of  fallback  in  each  ring  is 

then  obtained  by  dividing  the  fallback  mass  by  fpQ,  and  fallback 

volume  as  a function  of  range  is  found  by  summing  over  the  rings. 

In  scaling  results  to  different  yields,  fallback  must  be 
given  special  treatment.  The  position  and  time  of  departure, 
mass,  speed,  and  direction  follow  the  normal  cube-root-of-the- 
yield-ratio  scaling.  Since  the  ballistic  trajectories  are  deter- 
mined by  gravity  (independent  of  yield) , the  distance  and  time  of 
flight  do  not  scale,  but  are  constant  for  any  particle  (speed 
and  direction  of  exit  are  also  independent  of  yield ' . 

D3.4  Calculation  of  Bulking 

Bulking  of  material  outside  (below)  the  true  crater  re- 
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duces  crater  volume.  To  determine  the  amount  o£  bulking  that 
occurred  above  some  reference  depth,  D,  the  volume  of  material 
remaining  In  a given  ring,  after  ejecta  formation,  Is  determined 
by  subtracting  the  ejecta  volume  mej/pQ  from  the  ring  volume  D*A. 
The  volume  of  material  not  ejected  but  highly  cracked  (and 
bulked)  Is  then  found  from 

(D-D  ,)A  ; D>D  . 

ej  ej 


The  volume  of  the  material  after  bulking 


where  D . *A  = m ,/p  . The  volume  of  the  material  after  bulking 

is  then  V.  * gV,  , where  g is  the  ratio  of  the  initial  and  post 
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bulking  densities  pq/ . Studies  of  large  scale  craters  indi- 
cate that  g*1.05  is  appropriate. 

Two  different  values  of  the  reference  depth  were  considered: 
the  Depth  of  Burial  (DOB)  of  the  device,  and  a representative 
depth  of  cracking  (DOCRK)  for  a region  extending  slightly  beyond 
one  crater  radius.  The  first  is  considered  to  give  a more  rea- 
sonable estimate  of  the  effect  of  bulking.  The  second  overesti- 
mates the  effect  of  bulking  and  underestimates  crater  volumes. 

D4.  Crater  Dimensions 

Following  the  procedures  described  in  Section  D3,  contri- 
butions to  crater  volume  from  four  separate  mechanisms  were  com- 
bined in  the  following  ways  to  obtain  five  estimates  of  crater 
volume : 

1.  Ejecta 

2.  Ejecta  + plastic  flow 

3.  Ejecta  + plastic  flow  - fallback 

4.  Ejecta  + plastic  flow  - fallback  - bulking  above  DOB 

5.  Ejecta  + plastic  flow  - fallback  - bulking  above  DOCRK 


Case  2 represents  the  best  estimate  (in  this  set)  of  the  true 
crater  (if  it  had  the  same  shape  as  the  apparent  crater) , while 
Case  4 is  considered  the  best  estimate  of  the  apparent  crater. 
Plots  of  these  estimated  crater  volumes  as  functions  of  range 
are  given  for  each  of  the  9 basic  cases  in  Figures  D3  through. Dll 
Crater  radii  were  obtained  from  the  intersections  of  these  5 
curves  with  the  curve  for  the  volume  of  the  standard  inverted 
truncated  cone  crater  as  shown  typically  in  Figure  D2.  DOB, 
crater  radius,  scaled  DOB,  and  scaled  crater  radius  are  given 
for  the  five  estimates  of  crater  volume  in  Tables  D1  through  D5, 
respectively,  for  the  nine  problems  in  which  ejecta  cratering 
occurred  and  for  the  additional  cases  obtained  from  them  by  scal- 
ing. The  three  remaining  problems,  one  in  dry  sandstone  and  two 
in  wet  sandstone,  showed  that  no  ejecta  cratering  occurs  in  these 
media  for  scaled  DOB's  greater  than  or  equal  to  57.9m/ktA/  * . 


Scaled  crater  radius  data  from  Table  D4  (Case  4,  the  best 

estimate  of  apparent  crater  size)  has  been  plotted  as  a function 

of  scaled  DOB  in  Figure  D12.  Extrapolation  of  the  curve  for  NTS 

granite  shows  that  in  that  medium  no  ejecta  cratering  should  oc- 
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cur  for  scaled  DOB's  greater  than  about  80m/kt  ’ . 


Crater  Dimensions  for  a Revised  Method  of  Calculation 


Nine  calculations  of  ground  motion  from  tamped  150  KT  ex- 
plosions were  run  using  a revised  method  of  calculation  (Appendix 
C).  Four  of  these  were  recalculations  of  problems  which  produced 
ejecta  craters  in  the  earlier  set.  The  revised  calculation  pro- 
cedure affected  crater  dimensions  Scaled  crater  radii  from  re- 
calculation of  the  three  bursts  in  NTS  granite  showed  an  average 
increase  of  7.05%.  The  scaled  crater  radii  for  the  four  calcu- 
lations are  given  in  Table  D6.  Plots  of  crater  volumes  vs.  radius 
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Figure  D3 


Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
159.4m  in  NTS  Granite. 


nts  GftftNire 

•mm  w fULtmcn  it  o.im  « mm 

HUM  If  OUUtfO  HTCJtlOL  It  1.060  / 
tint  INTI  « Ito.t 
Hi  mi  • to?.* 

tOCf*  INI  s 303.0 
■ m ouuiihc 

• tuuiiMo  ftm  o oo 

a ouuiihc  rwn  oocot 

♦ cjrcva  oho  rvoouc 

K CJKCTO  ONLY 
O YOUHCOTtO  COHC 


Figure  D4. 


Comparison  of  Estimates  of  Cr  iter  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformatlon9  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
207.2m  in  NTS  Granite. 


.#  Figure  D6.  Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 

binations of  the  Computed  Contributions  from  Ejecta  Production, 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
, Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 

159.4m  in  Dry  Sandstone. 
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Figure  D7.  Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production, 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
207.2m  in  Dry  Sandstone. 
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Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
207.2m  in  Weak  Dry  Sandstone. 
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Figure  D9.  Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformation,  Fallback,  ar.d  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
53.1m  in  Wet  Sandstone. 


Figure  DIO.  Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production, 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
159. Am  in  Wet  Sandstone. 
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Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
207.2m  in  Wet  Sandstone. 


CRATER  RADII  ESTIMATED  FROM  COMPUTED  EJECTA  VOLUME 
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Figure  D12.  Scaled  Crater  Radius  as  a Function  of  Scaled  Depth  of  Burial 
for  Nine  Calculations  Cube  Root  Scaled  to  Yields  of  37.5,  150 
and  600KT  (Ejecta,  Plastic  Flow,  Fallback,  and  Bulking  above 
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based  on  the  five  estimates  of  crater  volume  are  given  in  Figures 
D13  through  D16  for  these  four  problems. 


Figure  D13.  Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
159.5m  in  NTS  Granite  (Revised). 
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Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters:  150KT  Burst  at  a Depth  of 
207.2m  in  NTS  Granite  (Revised). 


Figure  D14 


V/2PI  cn«*3)  Bio*4 

fl.H  -(•••  -Itl  IV. M tow  to. VI  40.00  10. 00  10  (1  70-00  09-00  00.00  100-00  110-00  ItO.OO  190.00  140-00 


MTS  GRANITE 

OCtftlTY  Of  r*LL6*C*  It  0.700  • KM61 

vetunc  or  ouuico  muoum.  it  i .os r / rhoi 

71*10  (NT)  e 150.0 

006  on  « *55.0 

00CJW  fin  c 460  *0 

O M0  OULMMO 

© OULU  t NO  IWfl  060 

a tuuiiNO  rxen  oocrx 
¥ tjccra  rno  ruisnc 

X CjeCTN  6NI.Y 
♦ TRUNCOTtO  C6H£ 


Figure  D15. 
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Figure  D16. 
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Comparison  of  Estimates  of  Crater  Volume  Based  on  Various  Com- 
binations of  the  Computed  Contributions  from  Ejecta  Production, 
Plastic  Deformation,  Fallback,  and  Bulking  with  the  Volume  of 
Standard  Truncated  Cone  Craters;  150KT  Burst  at  a Depth  of 
207.2m  in  Wet  Sandstone  (Revised). 
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APPENDIX  E 


f 


EFFECTS  OF  SURFACE  TOPOGRAPHY  ON  BODY-WAVE  MAGNITUDE 

Primary  calculations  in  this  investigation  addressed 
near- field  response  and  teleseismic  signals  from  buried  ex- 
plosions in  the  simplest  models  of  the  source  region  — homo- 
geneous half  spaces.  The  ground  surface  was  represented  but 
as  a simple,  featureless,  infinite  plane.  Computed  ground  mo- 
tions for  bursts  in  such  models  showed  the  major  effects  of 
the  ground  surface  — the  generation  of  reflected  signals,  inter- 
ference of  direct  and  reflected  waves,  and  the  development  of 
tensile  cracking  when  the  reflected  signals  were  of  sufficient 
strength.  The  latter  effect  accounted  for  significant  changes 
in  body- wave  magnitude  m^. 

Three  dimensional  surface  features  of  real  sites  — mou- 
tains , valleys,  canyons,  etc.  — within  ranges  comparable  to 
depth  of  burial  modify  the  motion  by  focusing  or  dispersing 
the  reflected  signals.  It  is  within  the  realm  of  possibility 
that  local  topographic  features  could  mask  or  significantly  mag- 
nify signals  from  buried  explosions.  In  the  parametric  investi- 
gation reported  here  the  effects  of  simple  topographic  features 
on  body-wave  magnitudes  and  m^  were  examined  for  bursts  in  a 
single  material. 

El.  Parametric  Study  of  Simple  Surface  Features 

To  focus  attention  on  the  influence  of  topography,  a 
number  of  simplifications  were  made.  Material  behavior  was  sim- 
plified: tensile  cracking,  shear  failure,  and  plastic  flow  were 

no  longer  allowed.  Linear  behavior  was  imposed  allowing  super- 
position and  simple  amplitude  scaling  to  be  used.  Topography 
was  simplified  to  permit  treatment  of  the  problems  in  2-dimen- 
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sional,  axial  symmetry.  The  resulting  site  models  consisted 
of  the  basic  homogeneous  half  space  modified  by  conical  fea- 
tures extending  out  to  a radius  of  1000m  from  ground  zero, 
where  they  meet  the  plane  surface  of  the  half  space  (Figure  E-l) 


Four  features  were  considered: 


20  Mountain  (Problem  4751.0) 


10  Mountain  (Problem  4753.0) 


Flat  Earth  (0  Mountain) (Problem  4750.3),  and 


10  Valley  (Problem  4752.0) 


Details  of  the  four  problems  are  described  below. 


El.l  Computational  Meshes 


The  finite  difference  meshes  consisted  of:  (i)  a stan- 
dard source  region  (Source),  (ii)  a region  of  slightly  skewed 
cells  (Region  I)  that  varies  from  problem  to  problem  to  account 
for  specific  surface  features,  and  (iii)  a basic  array  of  square 
cells  (Region  II).  These  are  shown  schematically  in  Figure  E-l. 
Details  of  Region  I and  the  Source  are  shown  in  Figure  E-2  for 
each  of  the  four  cases.  The  burst  points,  designated  by  stars, 
are  surrounded  by  five  points  designated  by  circles  and  squares 
marking  the  outer  corners  of  the  cells  in  the  standard  source 
region.  The  burst  point  (star)  is  always  253.0m  from  the  near- 
est point  on  the  surface:  the  circles  are  63.25m  away  vertically 
or  horizontally,  and  the  squares  are  63.25-  \7m  away  along  lines 
at  ±45°.  These  six  points  define  2 cells  63.25m  square.  The 
cells  in  Region  I are  then  obtained  by  expanding  (or  contracting) 
vertical  grid  point  spacing  between  the  source  cells  and  the  up- 
per and  lower  boundaries  along  both  the  axis  of  symmetry  and  the 
adjacent  line  of  grid  points.  Straight  lines  drawn  from  the 
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Figure  E2.  Source  and  Region  1 Mesb  for  Calculating  Ground  Response 
Affected  by  Local  Topographic  Features 
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column  of  cells  along  the  axis  to  the  corresponding  cells  at 
the  vertical  interface  between  Regions  I and  II  and  vertical 
lines  spaced  at  uniform  63.25m  intervals  define  the  remaining 
cells.  As  the  depth  approaches  2024m,  the  cells  smoothly  ap- 
proach the  uniform  63.25m  square  cell  configuration  of  Region  II. 

El. 2 Material  Model 

The  site  material  for  these  four  calculations  was  NTS 
granite.  In  this  series  of  problems  neither  tensile  cracking 
nor  shear  failure  was  allowed.  A linearized  version  of  the 
inelastic  model  of  NTS  granite  (Appendix  A)  used  for  the  matrix 
of  primary  problems  was  adopted  here.  In  the  linearized  version 
the  longitudinal  wavespeed,  C , is  4403  m/sec,  the  shear  wave- 

speed  is  2542  m/sec,  and  the  reference  density,  p , is  2.661 

3 ° 

gm/cm  . The  basic  63.25m  square  cell  provides  good  resolution 

of  signals  up  to  frequencies  of  10  Hz  in  this  material.  At  this 

frequency  P-  and  S-wave  wavelengths  are  about  7 cell  widths  and 

4 cell  widths,  respectively. 

El. 3 Source 

The  four  problems  in  the  parametric  study  were  driven 
by  velocity  histories  applied  at  the  five  source  points  (circles 
and  squares)  shown  in  Figure  E-2.  The  source  velocity  pulses 
were  defined  from  the  output  of  Problem  13,  spherically  symmetric 
motion  for  a burst  at  infinite  depth  in  NTS  granite.  The  reduced 
velocity  potential  (RVP)  from  this  problem  was  Fourier  trans- 
formed, low  pass  filtered  in  the  frequency  domain  to  remove  high 
frequency  content  (v>10Hz)  not  well  resolved  by  the  mesh,  and 
transformed  back  to  the  time  domain.  The  filtered  RVP  was  then 
used  to  obtain  radial  particle  velocities  for  ranges  of  63.25m 
(circles)  and  63.25-^2  (squares). 


Initial  attempts  using  these  source  velocity  histories 
showed  that  large  strains  developed  in  the  near  source  region. 

To  avoid  geometric  nonlinearities  associated  with  such  large 
deformations,  and  thus  allowing  invocation  of  linearity  and 
the  principle  of  superposition,  amplitudes  of  the  source  veloc- 
ity histories  were  reduced  by  a factor  of  100.  To  test  the 
validity  of  superposition,  the  flat  earth  problem  was  run  with 
(Problem  4750.5)  and  without  (Problem  4750.3)  an  initial  stable 
static  gravitational  field.  Dynamic  responses  showed  only 
minor  differences  in  the  near  field:  body-wave  magnitudes  for 
the  two  problems  derived  from  observations  at  angles  of 
0=70,75,80,  and  85°  showed  RMS  differences  of  only  0.005.  The 
variation  with  observation  angle  showed  an  RMS  value  of  0.064. 
Based  on  these  results,  the  four  problems  in  the  parametric 
investigation  of  the  effects  of  topography  were  run  without 
gravity. 

E2.  Results  of  the  Calculations 

Velocity  histories  were  obtained  at  1°  intervals  along 
a circular  arc  of  radius  3000m  centered  on  the  axis  of  symmetry 
at  the  elevation  of  the  planar  ground  surface.  Data  from  sta- 
tions at  0=70,75,80  and  85°  were  processed  to  obtain  seismometer 
displacement  histories  and  body-wave  magnitudes.  Body-wave 
magnitudes  m^  are  shown  for  each  observation  angle  for  each 
problem  in  Table  E-l.  Also  shown  are  the  average  magnitudes 
for  each  observation  angle  (averaged  over  the  four  topographic 
features)  and  the  average  magnitude  for  each  topographic  feature 
(obtained  from  the  average  maximum  seismometer  peak-to-peak  dis- 
placement A corrected  for  the  average  period  T over  the  four 
observation  angles).  Magnitude  variation  with  observation  angle 
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is  comparable  to  the  variation,  at  fixed  observation  angle, 
from  one  topographic  feature  to  another.  Comparison  of  the 
"average"  magnitudes  for  each  problem  shows  a magnitude  vari- 
ation of  only  0.068  over  the  four  topographies:  A difference 
of  0.235  occurs  between  the  20°  Mountain  and  10°  Valley  for 
the  75°  observation  angle;  a comparable  difference  of  0.186 
occurs  for  the  Flat  Surface  problem  for  observation  angles  of 
70  and  75°. 

The  overall  impression  from  these  results  is  that  there 
is  no  systematic  variation  of  m^  with  topography.  Examination 
of  values  of  m from  the  same  four  problems  and  their  averages 
by  viewing  angle  or  topographic  feature,  Table  E2,  shows  a 
strong  systematic  variation  with  topography.  For  all  four  ob- 
servation angles,  the  20°  Mountain  yields  the  lowest  value  of 
m . At  each  observation  angle  m increases  smoothly  as  the 
topography  changes  from  mountain  to  valley.  Using  the  "averages" 
over  observation  angle,  the  change  in  m in  going  from  the  20° 

cl 

Mountain  to  the  10°  Valley  is  0.458.  The  RMS  values  of  the 
variation  of  m with  observation  angle  for  fixed  topography 
and  with  topography  for  fixed  observation  angle  are  0.060  and 
0.400,  respectively,  for  m values  of  approximately  6. 

cl 

The  influence  of  topography  on  ground  response  is  more 
easily  seen  in  the  near-field  displacement  histories  and  tele- 
seismic  seismometer  displacement  histories.  Histories  from  the 
80°  observation  angle  show  typical  variations  with  topography. 

The  near-field  displacement  pulses,  Figures  E3-E6,  show  an  ini- 
tial positive  spike,  a negative  minimum,  and  a relatively  slow 
rise  to  a second,  slowly  decaying,  positive  peak.  The  relative 
amplitudes  of  the  peaks  and  the  shape  and  duration  of  the  nega- 
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Figure  E3.  Near-Field  Displacement  History  from  a Burst  Centered  253m 
Below  a 10°  Conical  Valley. 
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Figure  E5.  Near -Field  Displacement  History  from  a Burst  Centered  253m 
Below  a 10°  Conical  Mountain 
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Figure  E6 . Near-FieldoDisplacement  History  from  a Burst  Centered  253m 
Below  a 20  Conical  Mountain 
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tive  phase  indicate  that  the  amplitude  and  phasing  of  reflected 
signals  play  important  roles  in  determining  ground  response  at 
these  stations.  Without  interference  between  direct  and  re- 
flected waves,  the  initial  peak  should  be  decaying  as  r If 
the  differences  in  slant  range  due  to  variation  of  burst  eleva- 
tion relative  to  the  ground  plane  are  taken  into  account,  the 

-1  21  -2  38 

apparent  decay  rates  range  from  r to  r The  most 

rapid  decay  comes  from  comparing  the  two  mountain  topographies 
wherein  reflected  signals  should  be  focused  and  have  the 
strongest  influence.  The  weakest  decay  comes  from  comparing 
the  valley  and  flat  earth  problems,  the  case  of  diffusion  or 
weakest  reflected  signals. 

Also  very  evident  in  these  same  figures  is  the  marked 
increase  in  late  time  displacement  as  the  topography  changes 
from  valley  to  mountain.  At  t=2.1  seconds,  the  displacements 
are  2.1,  3.5,  4.8  and  6.0  cm  for  the  valley,  flat  earth,  10° 
mountain,  and  20°  mountain,  respectively. 

The  importance  of  interference  effects  is  even  more 
clearly  shown  in  the  plots  of  seismometer  displacement  histories 
in  Figures  E7-E10.  Shown  in  each  figure  are  the  peaks  AA'  and 
BB'  used  to  define  the  magnitudes  m^  and  m^,  respectively.  Even 
though  the  displacement  pulses  are  quite  similar  for  the  10°  Val- 
ley and  Flat  Earth  cases  (Figs.  E7  and  E8) , the  definition  of  the 
amplitude  for  m^  has  switched  from  one  feature  to  another.  The 
displacement  pulses  for  the  two  mountains  (Figs.  E9  and  E10)  are 
quite  different  from  those  for  the  valley  and  flat  earth  and  from 
each  other.  However,  the  feature  carrying  the  amplitude  used  in 
defining  m^,  though  different  from  that  used  for  the  valley  or 
flat  earth,  is  the  same  for  both  mountains. 
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Seismometer  Displacement  History  from  a Burst  Centered  253m 
Below  a 10  Conical  Valley 
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Figure  E8.  Seismometer  Displacement  History  from  a Burst  253m  Below  the 
Surface  of  a Plain  (0°  Mountain) 
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Figure  E9.  Seismometer  Displacement  History  from  a Burst  Centered  253m 
Below  a 10°  Conical  Mountain 
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In  all  four  cases  early  interference  phenomena  involve 
direct  (P)  and  reflected  (pP)  signals.  For  fixed  depth  of 
burial  the  delay  times  are  relatively  insensitive  to  topo- 
graphy so  that  the  variation  in  both  near-field  and  seismom- 
eter response  can  be  attributed  to  variations  in  the  strength 
of  reflected  signals.  This  is  controlled  by  the  focusing  or 
dispersive  effect  of  topography.  However,  altering  the  delay 
time  by  changing  depth  of  burial  or  the  medium  in  which  the 
burst  is  fired,  thus  changing  the  intrinsic  wavespeeds  as  well 
as  nonlinear  or  inelastic  behavior  affecting  signal  amplitude, 
may  significantly  change  the  relative  importance  of  topography. 


APPENDIX  F 
THE  SEDAN  EVENT 


F.l  Introduction 

To  establish  the  credibility  of  crater  dimensions 
determined  from  the  buried  burst  calculations  (Appendix  C) , an 
additional  calculation,  one  simulating  an  actual  cratering 
event,  was  carried  out.  The  Sedan  event  was  selected  from 
several  candidates  as  the  basis  of  the  demonstration.  It  was 
well  suited  for  this  role  since  1)  its  100-kt  yield  was  very 
nearly  that  of  the  standard  150-kt  yield  of  the  buried  burst 
problems,  2)  scaled  depth  of  burial  (50  m/kt  1/3.4)  fell  in  the 
middle  of  the  range  of  greatest  interest  (35  to  60  m/kt  1/3.4), 
the  site  consisted  of  a fairly  homogeneous  material,  4)  mate- 
rial properties  data  for  similar  materials  were  available,  5) 
crater  dimensions  were  well  documented,  and  6)  some  ground 
motion  data  were  available. 

The  demonstration  r*  quired  the  completion  of  several 
steps,  including:  1)  Description  of  the  test  event  (Section  F2), 

2)  characterization  of  material  behavior  (Sections  F3  and  F5), 

3)  development  of  a site  model  (Section  F4),  4)  Execution  of 
the  cratering  ground  motion  calculation  (Section  F6),  and  5) 
Evaluation  of  the  crater  size  prediction  (Section  F7).  Char- 
acterization of  material  behavior  and  development  of  the  site 
model  were  complicated  by  gaps  in  certain  of  the  data  and  by 
the  inhomogeneity  of  the  real  site.  These  "defects"  required 
some  departure  from  a pure  prediction  calculation,  i.e.,  one 
based  solely  on  preshot  material  properties  data  and  event 
geometry.  Equation  of  state  parameter  values  and  the  variation 
of  site  properties  with  depth  were  adjusted  iteratively  in  a 
series  of  preliminary  calculations  to  obtain  better  agreement 
with  observed  time  of  arrival  and  peak  velocity  data.  Finally, 
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the  shape  of  the  observed  Sedan  crater  was  used  in  place  of  the 
inverted  frustum  of  a cone  to  determine  crater  radius  from  the 
crater  volume  computed  from  ejecta,  fallback,  inelastic  volume 
change,  and  bulking  processes. 

F.2  The  Sedan  Event 

FI 

The  Sedan  event  was  a 100-kt  nuclear  cratering  experi- 
ment conducted  as  part  of  the  Plowshare  Program.  Ground  zero 
was  in  Area  10  of  the  Nevada  Test  Site  (NTS)  at  the  north  end 
of  Yucca  Valley,  near  the  sites  of  previous  cratering  events  of 
much  lower  yield  (Jangle  U,  Teapot  ESS,  Stagecoach,  and  Scooter). 

The  Sedan  device  was  placed  at  the  bottom  of  a 0.92-meter  (m) 
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diameter  cased  hole  at  a depth  of  193.5  m (50  m/kt  ).  At 
ground  zero,  the  alluvial  material  filling  the  valley  reaches 
its  maximum  depth  of  about  430  m.  A major  north-south  fault 
(Yucca  Fault)  is  located  approximately  1300  m west  of  ground 
zero.  West  of  the  fault  the  thickness  of  the  alluvium  is 
estimated  to  be  about  300  m.  East  of  ground  zero  the  depth  of 
alluvium  decreases  gradually  to  zero.  Variations  in  composi- 
tion and  cementation  have  been  observed  throughout  the  valley. 
Underlying  the  alluvium  is  a volcanic  tuff  with  physical  pro- 
perties much  like  those  of  the  alluvium.  Detonation  occurred 
at  10:00  a.m.  2 July  1962  producing  a nearly  circular  crater  of 
parabolic  cross  section  as  shown  in  Fig.  F2.1.  The  apparent 

crater  had  a radius  of  185  m,  a depth  of  98.5  m,  and  a volume 
6 3 

of  5.07x10  m . For  this  demonstration  its  shape  can  be 

3 x 

described  by  the  relations  V * 0.8007  R (or  R = 1.0769  Va)  and 
1 a a ' a a 

D = 0.5324  R . 
a a 
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Contour  map  showing  nearly  circular  symetry 


Figure  F2.1  The  Sedan  Crater 
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Figure  F2.1  The  Sedan  Crater  (concluded) 


F.3  Material  Properties 

Since  the  objective  of  the  Sedan  calculation  reported 
here  was  the  prediction  of  crater  size,  certain  liberties  were 
taken  in  preparing  the  site  model  and  the  equation  of  state  for 
the  alluvial  material  filling  the  valley  at  the  Sedan  site. 

Rather  than  follow  the  stringent  requirement  applied  to  ground 
motion  predictions  (which  prevents  use  of  any  data  from  the  test 
event  itself),  advantage  was  taken  of  the  experimentally  observed 
precratering  ground  motion  data.  Time  of  arrival  and  jump-off 
velocity  measured  near  ground  zero  were  used  to  evaluate  the 
adequacy  of  several  versions  of  the  alluvium  equation  of  state 
and,  later,  the  description  of  density  and  wavespeed  variations 
with  depth.  The  basic  form  of  the  equation  of  state,  described 
in  Appendix  A,  was  used  in  a series  of  calculations  of  one- 
dimensional, spherically  symmetric  ground  motion  to  iteratively 
obtain  an  optimal  description  of  Yucca  Valley  alluvium.  Equa- 
tion of  state  parameter  values  were  adjusted  until  satisfactory 
agreement  in  both  time  of  arrival  and  peak  particle  velocity 
were  obtained.  Then,  the  variation  of  density  and  wavespeed 
with  depth  were  varied  to  obtain  the  site  model  described  in 
Section  F.4.  Generalization  of  the  alluvium  equation  of  state 
developed  here  to  describe  alluvium  at  other  initial  densities 
or  wavespeeds  is  described  in  Section  F.5. 

F.3.1  Material  Shear  Behavior 

To  describe  the  shear  behavior  of  the  material  two 
quantities  are  needed:  The  shear  modulus  G,  and  the  yield 
surface  Y.  For  the  Sedan  calculation  the  alluvium  shear  modulus 
is  defined  by 
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G = MIN[G  , Q K U"  )]  [1-  (e/em)4] 
m u max  ra 
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where  G is  the  maximum  allowable  value  of  the  shear  modulus,  Q 
m 

is  a function  of  Poisson's  ratio  which  relates  the  shear  and 
bulk  moduli,  and  K (^  ' ) is  the  bulk  modulus  on  the  unloading 
hydrostat  at  the  compression  where  it  intersects  the  loading 
hydrostat.  The  quantity  em  is  the  energy  at  which  the  material 
behaves  hydrostatically  (i.e.,  it  melts  or  vaporizes).  The 
values  of  these  parameters  are  given  in  Table  FI.  Poisson's 
ratio,  based  on  the  unload-reload  curve,  is  0.466  and  the  max- 
imum value  of  the  shear  modulus  is  only  6.27  kbars. 

The  yield  surface,  Y,  is  a function  of  the  mean  stress, 
P,  and  the  specific  internal  energy  of  the  material,  e.  The 
yield  surface  is  given  by 

Y = Y(P)[1  - (e/em)4].  (F2) 

The  function  Y(P)  is  a standard  form  used  at  ATI  for  describing 
a generalized  Mohr-Coulomb-von  Mises  yield  surface. 


Y + Y'P 
c 


P * 0 


*(P)  - <YC  + (V^Kl-U-zn  p^  > P * o 


vm 


P 2 P 


vm 
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where 


TABLE  FI 

Material  Properties  Data  for  Sedan  Alluvium 


(bars)  -0.100000  b 0.6000 


and 


P 

vm 


YJ 
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For  samples  taken  from  the  Sedan  site  it  was  found  that  the 
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cohesion  of  the  material,  Y , was  about  2\  bars  . Data  on 

Merlin  alluvium  (from  NTS  Area  3,  several  kilometers  south  of 

the  Sedan  site  but  still  within  Yucca  Valley)  indicated  a von 

F3  i 

Mises  limit,  Y , of  180  bars  and  a slope,  Y , of  about  0.5. 

’ vm’  r * * 

The  pressure,  P , at  which  the  function  Y(P)  becomes  constant 
is  1420.0  bars. 

F.3.2  Mean  Stress  due  to  Compaction  in  Cold  Solid 

Since  permanent  loss  of  initial  porosity  must  be  accounted 
for  here,  the  equations  for  the  mean  stress  are  of  the  same  stand- 
ard form  as  the  equations  for  dry  sandstone  (Appendix  A).  The 
mean  stress  for  the  solid,  P , is  a function  of  the  excess 
compression  of  the  material,  ^ = p/pq  - 1,  and  the  maximum  value 

the  excess  compression  has  previously  attained,  p,  . If  n > p,  , 

max  max 

and  if  u s u , it  is  on  an 
max  ’ 

value  of  ll  . For  load- 
max 


n * ne 

^e  < u < uc  ’ 1 s s N 

P.  2 i±c  (F6) 


the  material  is  on  the  loading  cu 
unload-reload  curve  determined  by 
ing  the  mean  stress  is  given  by 

|Ki  “ 

ps=  y„D 

(F(u,tl*,VPp,Kp,Kbl> 
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r 


u. 


where  Sp(^,i)  are  cubic  splines  describing  the  P-4  curve  over 


N(=ll)  regions  of  the  excess  compression  interval  4e<4<4c>  and 


the  function  F is 


V * Km<“V  - (Km-Kp)>‘*{1-expC-(>‘-‘‘p)/'‘*^  +P, 


(F7) 


The  values  of  the  parameters  4*,  4 , P , K. , K , K are  listed 

’ p p ’ 1 * p * m 


in  Table  Fl  and  the  coefficients  for  the  cubic  splines  in 
Table  F2.  The  unloading  hydrostats  are  given  by 


[K.  ». 


max 


P = <P  (4,4  ) 

s \ u max' 


4 >4  and  4 < 4 

max  c c 


[F(4,4*,4p,Pp,Kp,Km)  4 * 4C  • 


(F8) 


The  function  P (4,4  ) represents  the  complicated 

u max 


intermediate  range  of  the  unioad-reload  hydrostat. 

The  general  form  of  the  unload-roload  hydrostat  for  the 


4 >4  and  4<4  , is 

max  e c ’ 


P = f P (4  ) 
p urvp 


where 


max 


11 


P (4  ) 

ur  'max' 


max 


S (4  , i) 

p pmax* 


P (x) 
ur  v ' 


[(K2  x + Kx)  x + Kt]  x 
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and 


pi  ( p,  ) = S (n  , i) 

^z'^max  p'^max5  ' 


The  relationship  between  the  functions  f and  pz  guarantees  P,p 

compatibility  at  the  load-unload  transition.  Fixing  either  f 

or  p,z  uniquely  defines  the  other.  For  instance,  at  a given 

value  of  u if  f is  specified  there  is  only  one  possible 
max  p 

value  of  p,  : conversely,  if  p,  is  defined,  the  value  of  f is 
z z p 

determined  uniquely.  It  is  numerically  easier  to  prescribe  p, 
and  solve  for  f^  since  the  expression  is  explicit.  Prescrib- 
ing fp  requires  an  iterative  procedure  since  the  expression  for 

p.  is  implicit.  The  function  S (u,  i)  is  a set  of  cubic 
z p ^max , ' 

splines  whose  coefficients  are  given  in  Table  F3.  The 
influence  of  changing  material  porosity  and  varying  seismic 
wavespeed  are  discussed  in  Section  F5. 


The  overburden  pressure  due  to  the  initial  gravitational 
field  is  added  to  Pg  to  obtain  the  mean  stress  for  solid  mate- 
rial. If  the  resulting  mean  stress  is  tensile  (negative)  and  is 

smaller  (more  negative)  than  a tensile  limit,  P , then  the 

’ ten  ’ 


mean  stress  is  set  to  P 


ten' 


That  is 


P 

s 


P 

s 


+ pb 


p 

s 


p 

s 


* p„ 

ten 


P 

s 


< P 

ten' 


When  P is  set  to  P , the  derivatives  of  the  pressure  with 
s ten  r 

respect  to  density  and  energy  are  set  to  zero  as  is  the  shear 


i*. 
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where 


modulus.  Since  this  material.  Yucca  Flat  alluvium,  is  not 
cemented,  the  tensile  cutoff  Pten  is  set  to  a very  low  value, 

-0.1  bars. 

F.3.3  The  Hydrostat  for  Sedan  Alluvium 

Although  no  pre-Sedan  hydrostat  data  are  available  for 

F 3 F4 

NTS  Area  10  alluvium,  some  hydrostat  data  are  available  ’ and 

F5 

data  for  Area  10  is  available  from  CIST-5  . Additional  data  are 

available  from  a "standard”  model  for  NTS  alluvium  used  by 

Shock  Hydrodynamics  (SH)  and  Physics  International  (PI)  in 

F6 

their  calculations  of  the  Johnie  Boy  Event  . An  early  equa- 

F7 

tion  of  state  developed  at  ATI  for  NTS  alluvium  was  also  used 
as  a basis  for  the  properties  of  Sedan  alluvium.  A hydrostat 
for  Sedan  alluvium  was  developed  from  a weighted  average  of 
these  data.  As  shown  in  Figure  F-3.1,  the  load  curve  has  a 
relatively  stiff  rise  to  about  50  bars.  This  is  followed  by  a 
region  of  high  compressibility  during  which  pore  closure  occurs. 
The  material  then  stiffens  as  the  pore  closure  is  completed  and 
the  compression  of  the  actual  grains  begins,  causing  the  fairly 
rapid  increase  in  pressure  with  continued  compression.  Per- 
manent compression  after  unloading  to  zero  pressure  from  rel- 
atively low  peak  pressures  (less  than  3.0  Kbars)  is  about  20-257. 

Since  much  of  the  data  on  which  the  Sedan  alluvium  model 
was  based  is  from  samples  from  other  sites,  the  model  was 
iteratively  tested  and  adjusted.  One-dimensional  (spherically 
symmetric)  calculations  of  the  precratering  motion  were  used  to 
monitor  "correctness"  of  the  model.  Parameter  values  were 
adjusted  and  the  calculation  repeated  until  differences  between 
calculated  and  observed  arrival  times  and  velocities  were  reduced 
to  a few  percent.  The  resulting  parameter  values  are  listed  in 


Table  FI.  The  cubic  spline  coefficients  S(ni)  for  the  load 
hydrostat,  and  S^(n  ) for  the  permanent  excess  compression 

lUaA 

are  given  in  Tables  F2  and  F3,  respectively. 

F.3.4  Energy  Dependence  of  the  Mean  Stress 

The  energy  dependence  of  the  mean  stress  is  separated 
into  two  distinct  parts.  The  first  accounts  for  the  effect  of 
heating;  the  second  describes  the  transformation  of  the  solid 
material  into  the  vaporized  state.  The  latter  is  complicated 
by  requirements  to  match  Hugoniot  data  on  shocked  alluvium  and 
the  behavior  of  the  vaporized  material  when  it  expands  to  large 
volumes . 

To  account  for  the  expansion  of  material  due  to  heating, 
the  excess  compression  used  in  the  equations  for  the  hydrostat 
for  cold  material  is  augmented  by  a thermal  expansion  compression 
term,  that  is 


+ e e 


where  0 is  equal  to  the  coefficient  of  thermal  expansion,  «, 

divided  by  the  specific  heat,  c . For  quartz,  the  basic  mineral 

^ F8  F9 

constituent  of  Sedan  alluvium,  the  values  of  a and  c are 

-6-1  v 
35.30x10  °C  and  0.188  cal/ (gm-°C)  respectively  which  give 


4.488x10 


gm/erg. 


The  high  energy  hydrostat  for  alluvium  is  of  the  well- 


known  Tillotson  form 


/ b e T) 

al+  a2  + — ~2 
\ e + e T\ 

' o 


exp(-»u/71)  . 
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The  original  set  of  constants  were  taken  from  an  earlier  equation 

F10 

of  state  for  alluvium  . These  constants  were  modified  to  give 

Fll 

better  agreement  with  Hugoniot  data  for  alluvium  and  to 

t 

increase  the  work  done  (and  the  strength  of  the  shock  driven 
outward)  by  the  vaporized  material  around  the  burst  point.  The 
changes  were  restricted  such  that  as  P-*-<®,  they  would  have  no 
effect  on  the  equation  of  state.  In  this  limit,  the  equation  of 
state  reduces  to  the  simple  form 

P = p(e  - ev)  (ax  + a2). 

To  satisfy  the  restriction,  the  quantity  (a^  + a2)  must  be  held 
constant  (=0.5).  The  parameter  is  equivalent  to  (y  - 1)  in 
an  ideal  gas  equation  of  state  where  volumes  are  several  times  that 
for  the  reference  state  of  the  material  (excess  compression  <0) . 
Changes  in  a^  and  a2  that  leave  their  sum  constant  do  not  affect 
the  Hugoniot  or  release  adiabats  for  positive  excess  compres- 
sions. Negative  excess  compression  states  (expansions)  are 
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affected.  Krieger's  data  on  Si02  were  used  to  obtain  a better 
value  of  a^  for  alluvium  vapor.  By  noting  that,  in  general,  the 
wavespeed  is  defined  by 


c 


2 


dP 

dp 


S 


and  that  for  an  ideal  gas  the  definition  reduces  to 

c2  = yP/p 
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r 


a simple  expression  for  y is  obtained 


dP/P 

d(ln  P) 

d(log  P) 

dp/p 

s ~ d(ln  p) 

s d(log  p) 

Replotting  Krieger's  data  in  the  form  P and  p vs  S (entropy) 
allowed  the  ratio  A(log  P)/A(log  P)|g  to  be  calculated  from  the 
values  of  log  P and  log  p at  given  S along  two  isotherms.  The 
values  of  log  P and  log  p on  any  isotherm  were  obtained  by 
logarithmic  interpolation.  Values  of  y from  9 points  were 
averaged  to  obtain  y = 1.18873  ± 0.01877:  all  values  fell 
within  the  range  1.16996  < y < 1.20750.  The  work  done  when 
using  the  value  y - 1 = 0.188773,  for  a^  with  a ^ - 0.5  - a^ 

= 0.31127  and  varying  the  parameter  b,  was  found  to  be  between 
58.77,  and  607,  larger  than  that  with  the  original  values  of  a^ 
and  The  work  was  greatest  when  b was  1.1.  This  value  of  b 

was  not  acceptable,  however,  since  the  resulting  Hugoiot  was 
stiffer  than  any  experimentally  observed.  Since  the  value  of  b 
does  not  significantly  affect  the  work  done,  the  value  which 
gave  the  best  fit  to  the  experimental  Hugoniot  data,  b = 0.6, 
was  used.  The  parameters  for  the  energy  dependent  part  of  the 
equation  of  state  are  given  in  Table  FI. 

F.4  Density  and  Wavespeed  Variation  with  Depth 

Since  the  alluvial  filling  of  a valley  is  a long  term, 
variable  process,  the  dependence  of  wavespeed  and  density  on 
depth  can  be  (and  usually  is)  more  complicated  than  the  simple 
variation  due  to  the  overburden  stress.  Pre-Sedan  seismic 
studies,  post-Sedan  borings,  and  additional  seismic  studies 
within  Yucca  Valley  indicate  a complex  variation  in  both 
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wavespeed  and  density.  Wavespeed  and  density  are  given  for 
several  depths;  details  of  the  variation  in  the  connecting 
regions  is  not  available. 

A piecewise  linear  fit  to  the  available  data  was  used  to 
define  the  initial  density  and  wavespeed  at  any  depth.  The 
available  data  did  not  justify  anything  more  complicated  than 
the  simplest  fit.  Since  the  points  of  the  transitions  are 
unknown,  within  the  limits  imposed  by  the  data,  their  location 
can  be  varied  to  correct  computed  arrival  times  and  jumpoff 
velocities  at  the  ground  surface.  The  initial  density  and  the 
wavespeed  ratio  c/c^  used  in  the  axisymmetric  calculation  of 
ground  motion  are  given  as  functions  of  depth  in  Table  F4  and 
plotted  in  Figure  F4.1.  c^  is  the  seismic  wavespeed  tor  the 
basic  loading  hydrostat  at  zero  excess  compression.  The  330  m 
depth  of  alluvium  in  the  site  model  adopted,  100  m less  than  the 
ground  zero  depth,  was  selected  as  representative  of  the 
region  (r  ~ 1 km)  surrounding  the  shot  point. 

F.5  Procedure  for  Varying  Porosity  and  Wavespeed  in  a 

Known  Constitutive  Description 

The  description  of  the  behavior  of  alluvium  given  in 
Section  F3  applies  to  material  of  specific  initial  density  and 
wavespeed.  The  site  description  and  site  model  (Sections  F2  and 
F4)  require  generalization  of  the  model  to  other  initial  den- 
sities and  wavespeeds.  The  material  properties  data  from  the 
Sedan  site  were  by  themselves  insufficient  to  define  the  model 
for  a single  initial  density/wavespeed  set.  Measurements  of 
the  hydrostatic  behavior  at  other  initial  densities  (air  void 
fractions)  are  not  available.  Generalization  of  the  alluvium 
model  to  the  range  of  observed  air  void  fractions  and  wavespeeds 
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TABLE  F4 


Depth  Dependent  Properties  in  Sedan  Alluvium 


Depth 

(m) 

°f 

(gm/cc) 

— 

C/Ci 

“of"  'f7'!-1 

Kof7Ki 

K * 
of 

(Kb) 

0 

1.700 

0.60 

-0.1052631579 

0.1733536156 

7.00048127 

33 

1.700 

0.60 

-0.1052631579 

0.1733536156 

7.00048127 

60 

1.835 

0.78 

-0.0342105263 

0.4841416195 

19.55092964 

90 

1.900 

0.78 

0.0000000000 

0.6084000000 

24.56881440 

104 

1.900 

0.78 

0.0000000000 

0.6084000000 

24.56881440 

166 

1.900 

1.00 

0.0000000000 

1.0000000000 

40.38266666 

245 

1.900 

1.00 

0.0000000000 

1.0000000000 

40.38266666 

280 

1.650 

1.00 

-0.1315789474 

0.3913664602 

15.80442130 

330* 

1.650 

1.00 

-0.1315789474 

0.3913664602 

15.80442130 

where  K^/K^  = (C/C 


) 2 i 

uof  I 


= 0.300796 

and 


K±  = 40.3826666  Kb 


Bedrock  is  at  330  meters  and  the  shot  point  is  at  195. 
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must  therefore  be  based  on  experience  with  other  media  and  some 
general  qualitative  constraints. 

The  following  are  examples  of  general  constraints:  (1) 
unload-reload  curves  for  a given  material  should  be  virtually 
independent  of  air-void-fraction,  (2)  curves  of  the  unload- 
reload  family  should  have  the  same  slope  as  the  loading  curve 
for  very  small  degrees  of  compaction  relative  to  the  overburden 
state,  (3)  the  loading  curve  should  be  tangent  to  the  limiting 
unload-reload  curve  where  the  two  curves  meet,  (4)  for 
undisturbed  material,  it  should  be  possible  to  set  the  slope  of 
the  loading  curve  to  any  reasonable  value  desired,  and  (5)  the 
loading  hydros  tat  should  have  no  more  than  a single  point  of 
inflection.  All  the  conditions  cited  can  be  met  by  simple 
transformations  of  the  P-  and  ^-axes  in  the  hydrostatic  stress- 
strain  plane.  The  transformation  described  below  meets  all  the 
stated  constraints.  It  therefore  provides  a means  of  cal- 
culating mean  stress  in  a medium  of  arbitrary  air-void-fraction, 
when  rules  for  calculating  the  mean  stress  in  the  same  medium 
are  known  for  one  specific  air-void-fraction. 

First,  various  quantities  that  appear  in  the  transforma- 
tion equations  must  be  defined: 

0 1 

a.  4 = elastic  part  of  the  cubical  dilatation  relative 

to  a single  reference  density,  pQ,  that  is 

independent  of  air-void-fraction 

i -Ael 

b.  = — r + Pe  = + Pe 

1 + 4 1 Z 

c.  = maxi-muin  value  of  reached  by  any  given 

material  element 
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d.  = value  of  ^ on  the  loading  curve  at  e = o and 

a mean  stress  equal  to  the  overburden  stress 


e . 


f . 


8- 


= value  of  where  the  limiting  unload-reload 
curve  meets  the  loading  curve. 

Pc  = value  of  P where  the  limiting  unload-reload 
curve  meets  the  loading  curve 

3P 

K ,.  = desired  value  of  r on  the  loading  curve  at  a 

of  *nf 


mean  stress  equal  to  the  overburden  stress. 

In  general,  the  basic  equation  for  defining  the  mean 
stress  P for  a given  general  state  is 


P (u  £ J U>  £ > I^q£  *^'of  ’ e ) = ^ (tA£  > U £ > t^of  ’^of  * e)  ?s 


,iAmax^ 


where  the  values  of  p,'  and  u are  functions  of  the  input 

^ max 

quantities  and  P (p,'  ,p,  ) represents  the  equations  defined  in 

s max 

Section  F3.2.  In  addition,  the  derivatives  of  P used  in  sound 
speed  calculations  are  with  respect  to  p,^  and  therefore  deriv- 
atives of  ^ ' and  cp  are  needed.  In  addition  to  the  quantities 
defined  above,  the  variables  K.  and  K , defined  as  the  values 
of  dP/dp,  at  the  overburden  stress  and  where  the  limiting  unload- 
reload  curve  meets  the  loading  curve,  respectively,  are  required. 

The  two  basic  quantities  needed  to  obtain  the  mean  stress 
for  new  void  fractions  and  wavespeeds  are  the  multiplicative 
factor,  cp,  and  the  transformation  to  an  equivalent  excess  compac- 
tion in  the  original  P-pt  space.  The  transformed  state  of  the 
material  depends  on  the  current  state  of  the  material,  p™j. 

If  nT  * or  p'  at  p,  then 
f c r e 
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Finally, 


dP 

9nf 


given  by  the  expansion 


, /p  iL  + -\  5P  54 

9p.f  \ dPm  / 


where 


1 

f 

9m.  £ 


1 


5P 

and  - — - is  the  slope  of  the  basic  hydrostat  in  the  state 
811 

(m'jM  v) . For  this  material  the  unload  slope  factor,  f , and 
max  p 

p.  (p.  ) are  obtained  by  iteration  to  insure  continuity  at  the 

z max 

point  of  intersection  of  the  unload-reload  curve  and  the  loading 
curve . 

F.5.  Constraints  on  Porosity  Model 

Experience  has  shown  that  constraints  are  needed  to 
insure  that  no  more  than  one  inflection  point  occurs  in  the  load 
hydrostat  for  the  porous  material.  These  constraints  are  given 
below , 


1)  • »c<‘‘of<l‘c 

11)  4(^c  + »of>Ki  > > ^c  + ,*of,Kl  ' °- 

Since  p,  = 0.300796,  K.  = 40,382.667  bars,  f = 2.93,  and  p„  = 1.90 
c i o 

for  the  basic  alluvium,  all  values  of  listed  in  Table  F4  are 


within  the  limits  set  by  condition  (i)  and  the  values  of 
found  from  the  table  satisfy  condition  (ii) . 

F.6  Calculation  of  the  Sedan  Event 

F.6.1  Preliminary  Calculations 

After  each  modification  of  the  equation  of  state  was  made 
and  tested,  a spherically  symmetric  (ID)  calculation  of  the 
motion  induced  by  a 100-kt  nuclear  device  was  made  to  determine 
if  the  resulting  motion  and  arrival  time  at  a point  193.5  m 
away  (the  ground  surface)  agreed  with  the  results  of  the  Sedan 
event.  As  part  of  this  series  of  calculations,  model  variations 
were  also  made  to  more  closely  simulate  the  axisymmetric  char- 
acter of  Sedan,  i.e.,  the  effect  of  gravity  and  the  variations 
of  density  and  wavespeed  with  depth  (radius)  were  introduced. 

At  least  one  aspect  of  the  model  was  changed  before  each  ID 
calculation  was  made.  After  about  10  calculations,  the  results 
were  in  good  agreement  with  the  experimental  data  and  it  was 
decided  that  no  further  refinement  was  needed.  Time  of  arrival 
was  within  2.5%  and  peak  particle  velocity  was  within  8.1%  of 
of  observed  values . 

Problem  execution  was  complicated  by  strong  variation  in 
density,  pressure,  and  energy  withii  the  vaporized  material 
filling  the  explosion  produced  ca^ ity.  Treating  the  vapor  as  a 
uniform  material  would  allow  significant  simplification  of  the 
calculation.  To  test  the  acceptability  of  this  change,  the 
spherically  symmetric  calculation  was  rerun  with  the  simpler 
model-uniform  vapor  driving  the  cavity.  The  results  of  this 
calculation  were  in  good  agreement  with  those  from  the  calcula- 
tion accounting  for  detailed  variations  of  pressure  within  the 
vaporized  material. 
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F.6.2  Axisymmetric  Calculation  of  Sedan 

Execution  of  the  2-D  axisymmetric  calculation  requires 
definition  of  a computational  grid  and  the  development  of  ini- 
tial conditions  for  material  in  the  cells  defined  by  the  grid. 

For  the  Sedan  calculation  initial  conditions  were  determined 
from  the  results  of  a 1-D  spherical  calculation  at  a time 
17.04  ms  after  detonation.  At  this  time,  the  initial  distur- 
bance has  not  yet  traveled  far  enough  to  encounter  material  of 
significantly  different  initial  density  or  wavespeed.  The 
cavity  has  expanded  to  the  point  where  no  further  material  is 
being  vaporized  and  the  uniform  pressure  model  of  the  vaporized 
material  can  be  used. 

The  initial  grid  for  the  Sedan  calculation  is  shown  in 
Fig.  6.1.  As  can  be  seen,  the  grid  in  the  half-plane  of  calcula- 
tion is  defined  symmetrically  about  the  horizontal  line  through 
the  burst  point.  The  region  of  symmetry  extends  vertically  and 
horizontally,  a distance  of  136.5  m from  the  burst  point.  Out- 
side the  rectangle  defined  by  the  two  136. 5-m  square  regions 
(above  and  below  the  burst  point),  the  grid  was  defined  by 
rectangular  cells.  Inside  each  square  the  grid  lines  transition 
from  a set  which  defines  the  surface  of  a quarter  circle  near 
the  burst  point  to  a set  of  vertical  and  horizontal  straight 
lines  at  the  outer  boundary.  The  initial  grid  can  be  divided 
into  the  three  regions  shown  in  Fig.  F6.2,  a blowup  of  the  grid 
in  the  vicinity  of  the  burst  point.  The  symmetry  of  this  por- 
tion of  the  grid  about  the  line  through  the  burst  point  is 
clearly  shown.  Region  I consists  of  material  not  yet  disturbed 
by  signals  from  the  burst  point.  The  initial  conditions  here 
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are  those  of  the  undisturbed  material.  In  Region  II  the 
material  has  been  disturbed  by  signals  propagating  from  the 
burst  point.  The  initial  conditions  for  these  zones  are 
obtained  by  combining  the  dynamic  field  variables  from  the  1-D 
calculation  with  the  field  variables  for  the  undisturbed  state. 
The  final  region,  Region  III,  corresponds  to  the  vapor-filled 
cavity.  The  uniform  pressure  assumed  to  act  in  this  region  is 
calculated  from  the  total  mass,  internal  energy  and  volume  for 
material  within  the  cavity. 

The  dynamic  field  variables  needed  in  Region  II,  are 
obtained  for  each  cell  by  evaluating  a weighted  average  of  the 
properties  in  the  spherically  symmetric  field  at  a time  of 
17.04  ms  over  three  range  intervals  defined  by  the  2-D  cell  as 
shown  in  Fig.  F6.3. 

As  the  problem  progresses,  the  spherical  symmetry  of  the 
initial  field  is  eroded  by  (i)  expansion  into  the  variable  earth 
and,  to  a lesser  extent,  (ii)  the  gravitational  field.  The 
symmetry  of  the  initial  state  is  shown  by  the  velocity  field  in 
Fig.  F6.4.  Some  ovalling  of  the  wave  front  is  evident  in  the 
velocity  field  for  130.3  ms  shown  in  Fig.  F6.5  (or  Fig.  F6.6). 

Throughout  the  early  stages  of  the  problem  the  cavity 
increases  in  size.  Near  the  cavity  wall  porous  alluvium  is 
compacted  to  about  4/5th  its  original  volume.  In  this  region 
near  the  cavity  wall,  the  radial  dimension  of  cells  experiences 
a large  decrease.  This  decrease  causes  a comparable  decrease  in 
timestep  and  the  calculation  becomes  inefficient.  To  allow  the 
problem  to  progress  at  a reasonable  pace,  zone  size  (and  the 
related  timestep)  are  increased  by  combining  the  thinnest  zones 
with  adjacent  zones.  Since  the  material  in  this  region  is  in  a 


DEPTH  in) 

-255-00  -235.00  -215-00  -195-00  -175.00  -155.00  -135-00  -1  15.00  -95,-0 

1 1 1 1 J I 1 J J 


RRNGE  ( li  ) 

•00  20-00  40.00  60-00  80-00  100.00  120-00 

1 1 1 1 l J 


Figure  F6.4  Velocity  Vector  Field  at  Start  of  2-Dimensional 
Axisymmetric  Sedan  Calculation 
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Figure  F6.5  Velocity  Vector  Field  at  130.297  ms:  Before 
First  Dezone 
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Figure  F6.6  Velocity  Vector  Field  at  130.297  ms:  Following 
First  Dezone 
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fairly  uniform  condition  with  no  appreciable  transient  motions, 
the  number  of  zones  can  be  reduced  even  further.  This  process 
of  combining  zones,  called  dezoning,  was  repeated  several  times 
during  the  calculation.  During  the  dezone  process,  momentum  and 
energy  are  conserved.  The  first  dezone  occurred  at  130.3  ms. 
Figs.  F6.5  through  F6.8  show  the  velocity  field  and  calcula- 
tional  grid  before  and  after  that  dezone.  Subsequent  dezones 
were  required  at  254.5,  286.2  and  351.6  ms. 

Changes  in  zone  size  and  shape  in  the  region  around  the 
cavity  can  cause  other  problems.  Zone  shape  was,  therefore, 
closely  monitored.  Severe  distortion,  which  may  lead  to 
instabilities  must  be  avoided.  By  utilizing  the  section  of 
AFTON  coding  which  allows  generalized  motion  of  cell  boundaries, 
minor  zone  shape  adjustments  were  made  over  many  cycles,  avoid- 
ing the  unwanted  zone  distortion  without  introducing  the  strong 
transients  associated  with  a single  sudden  zone  change. 

As  the  problem  progressed,  zones  were  added  radially  and 

vertically,  as  needed.  As  the  grid  grew  to  depths  greater  than 

330  m,  bedrock  material,  modelled  as  tuff,  was  added  at  the 

F13 

bottom.  The  Tuff  3 model,  developed  in  an  earlier  program  , 
was  used  to  describe  bed  rock.  As  can  be  seen  in  Fig.  F6.9,  the 
alluvium-bedrock  interface  causes  changes  in  the  appearance  of 
the  velocity  field.  By  450.7  ms,  (Fig.  F6.10),  the  motion  below 
the  interface  is  of  little  consequence  to  the  formation  of  the 
crater.  The  velocities  there  are  all  much  smaller  than  those  in 
the  crater  region.  No  sources  which  could  cause  strong  signals 
to  propagate  back  toward  the  crater  region  are  present,  below  the 
interface.  Since  only  weak  signals  will  appear  subsequently 
at  the  interface,  bedrock  can  be  replaced  by  a reflecting 
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Figure  F6.8  Computational  Grid  at  130.297  ms;  Following  First 
Dezone 
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boundary.  The  calculation  can,  therefore,  be  continued  without 
significantly  increasing  the  number  of  mesh  points  or  the  size 
of  grid  zones.  At  the  same  time  that  bedrock  was  replaced  by  a 
reflecting  boundary,  a radial  boundary  condition  was  imposed  at 
a distance  of  688.9  m from  the  axis  of  symmetry.  This  range, 
beyond  the  limit  of  signal  travel  at  the  time  the  boundaries 
were  introduced,  was  chosen  to  allow  use  of  the  grid  without 
further  alteration.  It  was  not  clear,  however,  which  type  of 
radial  boundary  condition  should  be  imposed.  Since  the  material 
would  not  behave  according  to  the  rules  of  linear  elasticity 
(permanent  compaction  takes  place  in  alluvium  at  very  low 
stresses),  an  elastic  transmitting  boundary  was  ruled  out.  Two 
options  remained:  (i)  a rigid  boundary  which  would  reflect  the 
entire  signal,  and  (il)  a floppy  boundary  which  would  simulate 
a free  surface  (i.e.,  the  motion  would  be  unopposed).  It  was 
decided  to  run  both  cases,  assuming  that  the  correct  result 
would  fall  somewhere  between  the  data  obtained  from  the  two 
calculations.  Both  problems  were  run  to  a time  of  2.666  sec, 
i.e.,  slightly  beyond  the  time  of  maximum  crater  radius. 

F.7  Calculated  Radius  of  the  Sedan  Crater 

As  the  two  problems  progressed,  the  radius  (volume)  of 
the  crater  was  determined  using  the  technique  described  in 
Appendix  C.  A slight  modification  was  made  in  the  definition 
of  the  theoretical  crater  shape  to  reflect  knowledge  of  the 
Sedan  crater.  The  relationship  between  crater  volume  and  radius 
was  evaluated  from  Sedan  crater  data  so  that  when  the  calculated 
crater  volume  was  that  observed,  the  deduced  crater  radius  would 
be  correct.  As  can  be  seen  in  Fig.  F7.1,  the  volume  of  the 
crater  continued  to  increase  until  a time  of  about  2.6  sec. 
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The  rigid  boundary  problem  crater, throughout  the  period  of  crater 

6 3 

formation,  was  about  0.2x10  m smaller  than  its  floppy  bound- 
ary counterpart.  The  actual  crater  radius  vs  time  curve  should 
fall  between  the  two  plotted  curves.  Tables  F5  and  F6  give 
radius,  depth,  and  volume  of  the  calculated  Sedan  craters  as 
functions  of  time  for  the  duration  of  the  calculation. 

Cratering  calculations  for  surface  bursts  typically 
predict  craters  whose  volumes  are  factors  of  30  to  50  low 
relative  to  experimentally  observed  values.  This  calculation 
resulted  in  a predicted  crater  volume  larger  than  the  experimen- 
tal results  by  only  56.870,  i.e.,  a factor  of  1.568  as  opposed  to 
1/30  to  1/50.  This  result  indicates  that  crater  radii  predicted 
by  the  technique  reported  in  Appendix  C should  be  in  close 
agreement  with  crater  radii  obtained  from  experiments  conducted 
in  the  materials  described  by  the  computational  models. 
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CALCULATED  CRATER  DIMENSIONS  AS  FUNCTIONS  OF  TIME 
RIGID  BOUNDARY  AT  A RANGE  OF  688.9  m 
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gravitational  effects,  indicate  a "reflection"  amplitude  smaller  than  that  cal- 
culated for  the  elastic  case.  More  important,  the  reflected  wave  is  relatively 
delayed  in  time,  so  that  transitions  between  positive  and  negative  reinforce- 
ment occur  at  shallower  depths  of  burial. 

A surface-wave  model  is  devc. loped,  based  on  Green's  function.  Many  prob- 
lems were  encountered  in  modifying  the  AFTON  source-data  program  to  provide 
information  that  was  accurate  for  long-period  displacements,  and  to  extrapolate 
calculations  well  beyond  the  reasonable  truncation  times  for  the  program. 

Preliminary  conclusions  are  made  concerning  the  need  for  inelastic  source 
calculations;  depth-of-burial  effects  on  signal  generation;  the  resulting  yield 
estimation;  possible  improved  yield  estimation  procedures;  and  topographic 
effects . 

Volume  I presents  summaries  of  the  body-wave  and  surface-wave  calculations 
to  date.  Volume  II  includes  additional  details  and  numerical  methods. 
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